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1. PROJECT SUMMARY 

This project is intended to research and develop new nonlinear methodologies for the control and 
stability analysis of high-performance, high angle-of-attack aircraft such as HARV (FI 8). The present 
report summarizes past research (reported in our Phase 1, 2, and 3 progress reports) and provides more 
details of final Phase 3 research. 

While research emphasis is on nonlinear control, other tasks such as associated model 
development, system identification, stability analysis, and simulation are performed in some detail as 
well. Table 1-1 provides an overview of various models that have been investigated for different 
purposes such as an approximate model reference for control adaptation, as well as another model for 
accurate rigid-body longitudinal motion. Only a very cursory analysis has been made relative to type 8 
(flexible body dynamics). Standard nonlinear longitudinal airframe dynamics (type 7) with the available 
modified F18 stability derivatives, thrust vectoring, actuator dynamics, and control constraints are utilized 
for simulated flight evaluation of derived controller performance in all cases studied here. 

Past research indicates that nonlinear control can be utilized effectively to control high-perfor- 
mance aircraft such as HARV (FI 8) for rapid maneuvers with large changes in angle of attack (a) — 
cases where classical linear feedback control (without gain scheduling) can yield poor performance or 
even instability. Even nonlinear feedback controllers that were generated in conjunction with a linear 
model reference (and without multiple regression terms) failed for certain high-a maneuvers but with 
added nonlinear reference terms they lead to successful control. More recent results, however, indicate 
that even the nonlinear feedback controller generated in conjunction with a higher-order (more delay 
terms) linear model-reference is quite effective for the cases considered in Section 6. On the other hand, 
it is shown that further performance gains may be realized with the addition of bilinear and higher 
nonlinear terms in the model reference for control adaptation as for previous studies. 

As a benchmark for judging controller performance, and also as a basis for future designs and 
for pilot-manipulated control, optimal control is studied in some detail, as summarized in Section 3. The 
derived algorithms for minimum -time control and for minimum quadratic state-error control both consider 
the full nonlinear dynamics with applicable control (magnitude and rate) constraints. The software has 
been generated in C language and can be applied virtually to any aircraft control or other plant. 

The time-optimal control from an initial trim state (at 15,000 ft altitude and mach 0.3) with 
afto) = 5° to «(%) = 60°, with thrust T(to) = 3,000 pounds and T(tf) = 18,000 pounds, takes place 
about 1.8 sec without any efforts to hold a(t) for t > tf . On the other hand, a quadratic-performance 
controller completes the maneuver in about 2.5 sec with a normal peak acceleration (n^ of about 2.5 g 
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Table 1-1. Aircraft Models 


Type 

Purpose 

Remarks/Limitatio ns 

1. Linear perturbations de- 

sired a, M, h 

Local control, check of nonlinear system, 
application of well developed linear con- 
trol methodologies 

Local stability 

Only valid for small maneuvers 
Special case of types 2-5 

2. Gain scheduled (non- 

linear function of a) 
from 1 

Gain-scheduled adaptive control based on 
well developed methodologies 

Simplified description of complex system 

Approximate stability 

May have stability problems with 
small number of reference states 
and/or large fast maneuvers 

3. Volterra seriest 

a) at reference states 

b) general case 

Nonlinear adaptive control via cross-corre- 
lation and/or £ priori dynamic structure 

Stability approximation 

Simplified dynamic description of complex 
system 

No n-orthogonal series approximation 

Sufficiency of 2 or 3 kernels 

Large computation time for adapta- 
tion 

4. Bilinear system 

a) continuous 

b) BARMA 

5. Polynomial time series 

Nonlinear adaptive control via model 
reference identification (NLMRAC) 

Stability approximation 

Simplified dynamic description 

Large computation time 

Bilinearizing controllers may be 
more practical than linearizing ones 

Polynomial approximation may be 
more accurate but more time con- 
suming than linear or bilinear ap- 
proximation 

6. Neural network 

Potential application to adaptive control 

Probably less accurate than 4 or 5 
for a given data set but accuracy 
may be more robust outside the 
available data set 

7. Nonlinear ordinary dif- 

ferential equation model 

Accurate approximation to fast large ma- 
neuvers for "final" design and simulation 

Stability 

Neglects flexible modes and other 
complications 

8. Nonlinear partial differ- 

ential equations model 

Allows the treatment of distributed exter- 
nal stresses and local pressures as well as 
internal stresses and deformations 

Computations are time consuming 

Overall model complexity 

Allows one to treat flutter and dis- 
tributed control strategies, etc. 


*ct is angle of attack, M mach number, and h altitude 
f Wiener series can be used for orthogonal representation 
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and a q max (maximum pitch rate) of about 60°/sec. The a response is quite smooth with less than 10% 
overshoot and a rise-time of about 1.8 sec (see Section 3.4). 

All the high-performance controllers studied seem to mimic the optimal control policy by 
operating in saturation for most of the trajectory. Switching times of the stabilator and the thrust vector 
motion is somewhat critical as generated by the appropriate feedback. 

The stringent requirement to control a and not allow drastic changes in pitch for t > % suggests 
the use of a terminal sliding-mode control which is studied for die whole maneuver in Section 7. While 
the results are preliminary, the method shows promise. 

The adaptive controls studied in Section 6 perform very well and only a little slower than 
optimum. It is expected that further investigation of practical nonlinear feedback control will lead to an 
even simpler and better closed loop system. Then effectiveness for trajectory following in general needs 
to be studied. 

The neural-net controller which is trained off-line on the simulated optimal policies perform quite 
effectively despite the preliminary nature of derivation. For example, such neural-net generated policy 
controls HARV (F18) from a a= 5° trim state to 70° in about 4.5 sec. It should be realized that, in a 
general sense, neural-net-based controllers typically are a form of nonlinear adaptive controllers with the 
neural-net an appropriate nonlinear model consisting of layers of adjustable parameters with their outputs 
summed and processed through nonlinear limiting functions such as sigmoidal ones. The other adaptive 
controllers here are based on polynomial (with linearity a special case) discrete approximations with 
varying degrees of lag terms. Project publications (in addition to the semiannual progress reports) are 
listed in Appendix A along with research contributors. 
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2. INTRODUCTION 

Nonlinear control of high-performance aircraft such as HARV (FI 8) has been demonstrated with 
accurate computer simulations to be effective by this project as well as others such as Ostroff [2.1] and 
Buffington, Sparks, and Banda [2.2]. The difference between our controllers and the latter two is that 
the latter generate the nonlinear feedback gains to account for certain (assumed known) nonlinear plant 
dynamics. Reference [2.1] utilizes numerous trim-state linearization studies to determine the required 
gains in conjunction with a flight-test proven PIF controller which is able to successfully control HARV 
from a(to) = 5° trim to about 60° in about 5 sec. Reference [2.2] is based on a linear H* design in 
conjunction with trim-state linearized dynamics and an appropriate nonlinear gain scheduled according 
to dynamic pressure variation. The latter study only considers a maximum change from a(to) = 10° to 
a = 20° in about 3 sec with a rise time of 1 sec, q,,,^ ~ 14°/sec, and ^ — 1.5 g. While neither of 
the latter two are nearly minimum -time maneuvers as demonstrated here, they probably represent the best 
controllers based primarily on linear design methodology in conjunction with somewhat ad-hoc nonlinear 
correction. 

2.1 References 

[2.1] A. J. Ostroff, "High Alpha Application of Variable-Gain Output Feedback Control," J. 
Guidance, Control & Dyn., Vol. 15, 491-497, 1992. 

[2.2] J. Buffington, A. Sparks, and S. Banda, "Full Envelope Robust Longitudinal Axis Design of a 
Flight Aircraft with Thrust Vectoring," Proc. IEEE American Control Conf., San Francisco, 
1993. 
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3. OPTIMAL CONTROL 
3.1. General Concept 

This section discusses research on optimal control of longitudinal aircraft dynamics. Since the 
objective of the project is to develop high performance control algorithms, it is natural to investigate the 
possibility of obtaining the highest possible performance — through optimal control theory. While 
optimal control may be useful in its own right, the emphasis here is in using the optimal control as a 
yardstick to judge our more practical designs and as a base for their synthesis such as by neural nets. 
Understanding how the formulation of the optimal control problem affects the behavior of optimally 
controlled aircraft is crucial to a successful application. In this work two most commonly used 
performance criteria are investigated: minimum time performance index and integral quadratic 
performance index. 

Even if the correspondence between the problem formulation and actual optimal performance is 
well understood, still practical application of optimal control theory sometimes may seem questionable. 
Most significantly, optimal control depends upon exact knowledge of the dynamics of the controlled plant. 
While this assumption may be quite realistic in some applications, parameter uncertainty arises for high 
performance aircraft due to stability derivative uncertainty, nonuniform (geometrically) fluid flow, etc. 
Also, for a nonlinear system, the actual calculation of the optimal controls usually involves an iterative 
numerical algorithm. Even in case of a relatively uncomplicated model, as longitudinal aircraft dynamics, 
on-line optimization may require computational capabilities of the order precluding any realistic practical 
applications. Furthermore, unlike for linear systems, optimal feedback synthesis for a nonlinear plant 
may be almost impossible. There are numerous optimization techniques available that solve open-loop 
optimal control problems efficiently. However, there are few methods to solve the closed-loop synthesis 
problem. Theoretically, the very existence of a unique optimal feedback is not even certain for some 
nonlinear systems [3.6]. These three points make one raise a question: If the available aircraft model 
is highly uncertain and feedback control policy is sought that must be implementable in real-time using 
an on-board computer is it at all logical to investigate optimal control? It is argued here that in spite of 
those doubts the study of optimal controls can answer important questions and be utilized in the design 
of feedback controllers. 

First of all, optimal trajectories provide the upper bound for the performance possible achieved 
by the investigated aircraft. Since the solution of the optimization problem is obtained with full and exact 
information about the plant dynamics, no other control algorithm can possibly result in better performance 
in terms of a particular criterion. Thus the optimal trajectories can serve as benchmark tests for feedback 
controllers developed using other design methodologies. With the limit of performance available, it is 
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easier to assess the degree of improvement of one control scheme over another. This alone seems to be 
an objective important enough to justify the investigations that will follow. 

Second, optimal trajectories demonstrate what type of control actions are required to achieve high 
performance. In particular, optimality is often achieved with extremal values of control signals. Current 
practice in aircraft control design is to use gain scheduling of linear controllers, and often the point is 
made of avoiding the saturation of actuators (e.g. [3.4]). This is motivated by the fact that saturation 
results in loss of linearity of the controller and invalidates the theory underlying the design. The issue 
emphasized here is that, since the aircraft dynamics during large maneuvers can hardly be described by 
a linear model, the struggle to maintain linearity of the controller is somewhat unnecessary. As it is 
demonstrated below, almost always maximum admissible control effort is required for both time-optimal 
control and integral-quadratic criterion minimization. Thus, to utilize existing control capabilities fully 
to their limits, it is necessary to use nonlinear control algorithms, which take into account limits of 
control signals and of their change rates. Investigations of optimal open-loop trajectories may also answer 
the questions which of those constraints are most critical for which types of maneuvers. This may 
provide some guidelines for the designers of the actual aircraft about possible ways of enhancing optimal 
aircraft performance. 

Third, open-loop generated optimal trajectories may be used for synthesis of feedback controllers. 
One such possibility, discussed in detail in Section 4, is to approximate optimal feedback mapping using 
a n umb er of optimal trajectories as learning data used for nonlinear function fitting. A particular 
imp lementation of such approximation proposed in this report is by artificial neural networks. Another 
way to utilize information about optimal trajectories is to use them as reference tracking signals in output 
feedback controllers. Formation of reference signal out of command signal is one of the major difficulties 
in model reference control. Common practice is to use linear reference models. As discussed above this 
may lead to unnecessarily conservative performance, particularly close to the target value of the output. 
A solution to this problem may be the utilization of optimal trajectory as the demanding, yet possible to 
track, reference signal. Difficulty with on-line generation of an optimal output trajectory may be 
circumvented by some kind of interpolation between a number of optimal trajectories, pre-calculated off- 
line for different command signals. Again, the artificial neural network approach is one possible 
implementation of this idea. 

To investigate the issues raised above a large number of optimal trajectories was calculated in the 
open-loop mode using numerical optimization techniques. The following sections discuss in detail optimal 
control problems solved, the actual algorithms used for solution, and the results of simulations. 
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3.2. Methodology 

The optimal control calculations were performed for the fourth-order nonlinear model of 
longitudinal dynamics augmented with dynamics of three actuators including their saturation. Two types 
of optimal control problems are solved, time-optimal problem and integral quadratic performance index 
minimization, with controlled output variable being either the angle of attack or pitch rate. The following 
subsections give precise problem statements and discuss the numerical methods used. 

3.2.1. Aircraft Model 

The model of the aircraft used is the full nonlinear longitudinal dynamics of a modified F 18 
(HARV), described in Appendix B. The general form of the models is: 

dx/dt = f(x,u) @.l) 

with state vector x = [a, q, 6 , v] T , control vector u = [5 h , T M , S V ] T , a — angle of attack, q — pitch 
rate, 6 — pitch angle, v — total velocity, S h — elevator angle, — thrust magnitude, and S v — angle 
of thrust vector. In practice, thrust magnitude is controlled typically by the pilot through the throttle, 
so generally it should be treated as a given, time-varying parameter, rather than a feedback control signal. 
Here, the objective is to obtain the performance limits for the aircraft and to manipulate T M optimally. 
Control signals are assumed to be subject to the following constraints: 

-24° < 6b <£ 10° 

3,000 lbs <£ T m < 18,000 lbs (3.2) 

-20° < 6 V < 20° 

Additionally, model (3.1) was extended with dynamics of the actuators for 5 h> t m> and 5 V . Dynamics 
of T m is assumed to be linear of the first order with time constant Is. Time constants for 5 h and <5 V are 
l/30s, but there are constraints for their rates of change: |dS h /dt| < 40°/S, | dSy/dt | < 40°/S. This 
results in nonlinear dynamics model for Sj, and 6 V . The saturation of rate of change is modelled in a 
smooth way to make the gradient calculations in optimization methods possible. The actuators were 
modelled as: 

dydt =g(-305 h+ 30u 5h ) 

dS v /dt = g(-306 v + 30u 5 J (3-3) 

dT M /dt = -T m + u Tm 
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with smooth saturation function g 


g(z) = 


40 -exp(-x+39) 
z 

-40 + exp(x+39) 


if x > 39 
if -39 < x < 39 
if x < -39 


Finally the obtained model is of the seventh order 


dx'/dt = f'(x',u') 


(3.4) 


(3.5) 


with the augmented state x' = [x T , u T ] T and with the control variables being the demand signals to the 
actuators u' = [u 5h , u T m , u 5y ] T . Constraints (3.2) are now defined for u', but because of the form of 
(3.3-3. 4) will be also satisfied for actual control signals u. 

For the purpose of study of optimal control, equilibria of model (3.5) were investigated. Since 
in equilibrium u = u', it is possible to omit three equations (3.3) and deal only with the fourth order 
system (3.1). Because thrust vectoring is usually applied only during maneuvers, 6 V in equilibrium is 
assumed to be zero. In such case the set of equilibria of (3.1) is parametrized by values of T M and <5 h 
— it is a two-dimensional manifold in four-dimensional state space. Among all equilibria there are trim 
conditions, e.g., specific points in which a = 6. For each thrust magnitude T M there is only one or two 
such trim conditions, and all other attainable equilibria of (3.1) correspond to ascending or descending 
flight. For higher T M it is possible to trim the aircraft at higher a — for maximal thrust of 18,000 lbs 
maximal trimmed a is close to 35°. For angles of attack or pitch angles larger than 35° there are 
equilibria of (3.1) but with non-zero climbing angle. 


3.2.2. Time-Optimal Control 

The control problem is as follows: given initial state x° and desired final state x f minimize 

J to = min{t f :x(t f ) = x f } (3 ' 6) 

with respect to control signal u(t), 0 < t < / . In what follows the plant considered is (3.5) even though 
primes are omitted to simplify notation. The minimal time problem is solved in hierarchical manner. 
The first step is the problem of reaching x f from x° in fixed time t*. To solve this an auxiliary 
performance criterion J aux is minimized with respect to control signals: 


3-4 


(3.7) 


Here p x are appropriate weights playing the role of scaling coefficients. Also, setting some of to zero 
allows us to specify the desired final set as a hyperplane rather than a point. This may be the case if 
values of only some of the state variables are actually of interest. If the desired state x^ (or, more 
generally the hyperplane specified by and p) is reachable from x° in time t* then the minimum of (3.7) 
is J^(t f ) = 0. Therefore the actual time optimal problem (3.6) can be reformulated as a search for: 

t°P' « min{t f : J^(t^ = oj <3 8) 


Minimization of (3.7) is performed using an iterative gradient scheme. The key element of the 
procedure is the calculation of the gradient of J aux with respect to control signals. This involves solving 
at each step a system of adjoint equations 

djj/dt = -(9f(x,u/dx) T ij (3 ’ 9) 

with final conditions 


i?i(t f ) = -2pi(xj(t f ) -xjf) 


(3.10) 


Then the gradient can be calculated according to: 

<aj aux /9u)(t) = ((df/duVvh 


(3.11) 


For derivation of these sensitivity calculations see [3.7]. The gradient calculated according to (3.11) is 
used for determination of change direction in which minimum search is performed. In this study a 
conjugate gradient method (see e.g., [3.2]) is used for this purpose. Simplified flow of information 
during each iteration of minimiz ation of (3.7) is depicted diagrammatically on Fig. 3.1. The most time- 
consuming element of each iteration is the directional search which consists of repeated simulation of the 
model with control signal u(t) + yv(t), where y is the search step. Implementation issues concerning 
directional search, control discretization and incorporating constraints for control values are discussed in 
Section 3.2.4. 
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new u(t) 

Figure 3.1. Flow of Information for Minimization of (3.7) 


With the method of minimiz ation of (3.7) available, determination of optimal time (3.8) is a 
simple one dimensional optimization problem. In this study it is solved using a bisection method which 
for our purposes proved to be effective enough. 

An alternative, and often more numerically efficient, algorithm for minimization of (3.7) is the 
switching time variation method [3.3]. It seeks the solution in form of bang-bang controls with the 
independent variables being the switching times. If optimal control is indeed of bang-bang type with 
finitely many switchings then the resulting optimization problem is usually low-dimensional and converges 
mnrh fester than optimization with respect to control values. If the time optimal control trajectory 
contains a singular arc (an interval with control values not on the constraints) the switching time variation 
method is still applicable, since for plants affine in control any control trajectory can be approximated 
arbitrarily closely with bang-bang control [3.1]. Technically (3.3) is not affine in control, but is strictly 
increasing and by trivial transformation may be made affine in a modified control signal with one-to-one 
correspondence to the actual control signal. Unfortunately, in a situation with a singular control arc, the 
method results in increasing number of switchings and quite slow convergence. On the other hand 
optimization with respect to control values, rather than switching times is more general and also suitable 
for integral quadratic performance criterion optimization. Therefore application of this method made it 
possible to use the same basic optimization module for different performance criteria. 

3.23. Optimal Integral Performance Criterion Control 

Control trajectories resulting from solution of the time-optimal problem often display large over- 
shoots which might be unacceptable from the designer’s point of view. Also they include no information 
as to how to keep the final state. In many time-optimal control applications this issue is not relevant, as, 
e.g., in case of missile guidance. If, however, the control task is to change the angle of attack or pitch 
angle rapidly it is expected that the angle should be kept at the same level for at least some finite time. 
The problem with terminal end control is most apparent if only some of p, are non-zero — i.e., if values 
of only some of the states are of interest. Then the time-optimal movement towards the target hyperplane 
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may result in such a final state for which no admissible control exists which would keep the state on the 
hyperplane. This difficulty can be overcome by using integral performance criterion: 



(3.12) 


A typical choice of the integrand is a quadratic error function 

f 0 (x,x f ,t) = 

If the final state x f is an equilibrium corresponding to an admissible value of control then intuitively one 
may obtain the steady state control by minimiz ation of (3.12-13) with ^ -* oo. Even if is not an 
equilibrium, still optimization with finite but large t f yields a control which keeps the state close to the 
desired value. If the final set is not a point but a hyperplane (some of Pj are zero), minimization of 
(3.12-13) may result in controls keeping particular state coordinates exactly at desired values. This is 
the case for the control problems discussed in this study. 

Criterion (3.12) is minimiz ed in an iterative scheme, whose main elements are: calculation of 
current value of the criterion, calculation of its gradient with respect to the control signal, determination 
of the search direction and directional minimiz ation. To find the gradient it is necessary to solve 
backwards in time the system of adjoint equations: 

dif/dt = -(3f(x,u)/dx) T i/ + (df 0 (x)/dx) T (3-14) 

with final conditions 




(3.15) 


Then the gradient of J with respect to control signal u is calculated according to formula (3.11). The rest 
of the optimization process is also similar to that discussed in Section 3.2.2 and depicted on Fig. 3.1. 


3.2.4. Implementation 

The actual numerical optimization of (3.7) and (3. 12) in this study was performed in discrete time 
— with control signal piecewise constant on intervals with length A (except possibly for the last interval, 
which may be shorter if / is not commensurate with A). Discrete time approximation allows for 
utilization of all the powerful machinery of finite dimensional optimization. The infinite dimensional 
control signal u(t), 0 < t < t f is replaced with finite dimensional vector u = [uj, u 2 , ...» j ] T 
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0 < k < K f , where K f is the smallest integer such that t f < K f A. Consequently, the gradient of the 
performance criterion with respect to control signal is approximated by a finite dimensional vector g 
whose elements are calculated according to: 

gk - *>W 3u k ■ |»* ,)a ((9f/8u) T ,)dt (3 - 16) 

Then the directional minimum search is performed in direction v determined by a conjugate gradient 
algorithm. Here a particular version of that algorithm is adopted from [3.5], which is robust with respect 
to inaccurate direction of search. This allows us to limit the number of calculations of the performance 
criterion (3.7) or (3.12) in each iteration. 

A pure conjugate gradient algorithm is suitable only for unconstrained problems — in this case 
with no limits on control values. With constrained controls the minimum in direction v may be infeasible. 
This poses no difficulties in the directional search — in fact it is easier to find a minimum with a 
constrained search interval. However the resulting optimal step cannot be used for determination of the 
next search direction if the minimum lies on one of the interval’s ends. This problem is solved as 
follows. After calculation of (3.11) the set of active constraints is calculated: the constraint is active if 
g k < 0 and u k = or if g k > 0 and u k = u,^. Then the whole minimization is done only with 
respect to controls with no active constraints — i.e., elements of g corresponding to active constraints 
are set to zero. If some of the active constraints become non-active or vice-versa, then the conjugate 
gradient algorithm is re-initialized with the steepest descent direction as the initial search direction. The 
same is done if in some iteration the search direction v becomes infeasible. As a result the conjugate 
gradient method works effectively as if no constraints were present and the formulae for determination 
of search direction are valid. Hitting or leaving any of the control constraints causes change of the space 
in which the minimizati on is performed. The method used may be viewed as a projection method ([3.2]) 
in which currently active constraints are treated as equality constraints and the others are ignored. Due 
to the extremely simple structure of constraints (3.2) the projection of the gradient amounts to putting 
some of the gradient elements to zero. Also, since the constraints are linear, they are satisfied for any 
control signal along the search line. 

The directional minimization implemented in this study uses a combination of two-point gradient- 
based parabolic fit with three-point non-gradient parabolic fit. If the calculated vertex of the fitted 
parabola results in controls exceeding the constraints, a heuristic bisection-based search for minimum in 
the interval is performed. The number of criterion value calculations (i.e. simulations of the model) in 
each directional search was limited at cost of accuracy. 
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The preceding discussion, for sake of clarity, was limited to the scalar control case. Extension 
to more control signals (such as the case of aircraft control) is straightforward. The control vector, with 
respect to which the minimization is performed is now u = [u{ 1 ^...,u{ M ^,U 2 ,--.,u^ M ' 1 \...,i 4 ^-iv, 
ujjfif, where M is the number of control signals. The gradient of the performance criterion is composed 
of M vectors g = [g (1)T ,...,g^ T ] T , each of which is calculated according to (3.16) with u being 
replaced by appropriate control u^. 

For the time-optimal problem the minimiz ation of (3.7) is terminated when J aux < e, where e 
represents acceptable deviation from desired values, or, if for given t* it is not possible to reach x f , when 
the norm of gradient of J becomes small enough. For the integral-quality-criterion problem, only that 
second stop criterion was used. 

3.3. Simulation 

For the aircraft model a series of numerical optimization experiments was performed for both 
time-optimal and quadratic integral criterion problems with controlled variable being either angle of attack 
a or pitch angle 6. The original objective for the time-optimal problem was to calculate the minimum- 
time transition of the plant’s state from one point (assuming it is an equilibrium of (3.1)) to some other 
specified point (preferably also an equilibrium). For the control of or (or 6) the procedure would be to 
find an equilibrium corresponding to a desired value of a (or 6 ) and then calculate the time-optimal 
trajectory. Unfortunately this turned out to be too difficult a problem with given control limits, and it 
was impossible to find appropriate controls moving the state to such a desired equilibrium in a reasonable 
time. Besides that, equilibria corresponding to high values of or are characterized by negative values of 
path angle, which is not exactly the desired state. Therefore a modified control problem is solved, in 
which the task is only to reach a given value of or (or 6) and q = 0. This corresponds to = p 4 = 0 
(or pj = p 4 = 0) in (3.10) and to the target set being a two-dimensional hyperplane rather than a point 
in R 4 . This turned out to be easily solvable; however, specifying final values of three, instead of two, 
state variables caused serious problems with convergence of the minimization process due to a physical 
lack of aircraft controllability with present control signal constraints (particularly thrust). 

Typical time-optimal trajectories are displayed on Figs. 3.2-3.6. The initial condition is the 
trimmed flight with T m = 3,000 lbs at 0.3 mach and 15,000 ft. The corresponding pitch angle is close 
to 5°. The target is q = 0 and 6 = 10°, 30°, 50°, 70°, so the simulated maneuver corresponds to 
aimin g the nose of the aircraft in the desired direction. The dashed lines on control plots correspond to 
command signals fed to actuators while the solid lines are the actuators’ outputs, i.e. the control signals 
proper. It may be observed that the optimal command signal fed to the actuators is not unique because 
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Figure 3.5. Time-optimal control to 8 = 70° and q = 0°/s 
Initial condition - trim corresponding to thrust 3,000 
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of the rate limitations on actuators’ outputs. Theoretically the nonlinear saturation function g in (3.3) is 
monotone increasing, so bang-bang co mm and signals are necessary to obtain the fastest changing control 
signals. However, from the practical point of view, g saturates so quickly, that any command signal 
differing from current actuator output by more than 2° results in the maximal rate of change. Results 
displayed on Figs. 3. 2-3.5 show that indeed time-optimal change of pitch angle requires bang-bang 
control of all three actuators. This substantiates previous intuitive remarks about unsuitability of linear 
methods for high-performance control. Optimality requires utilizing the existing control capabilities fully 
to their limits, and increasing the speed of actuators is one of the ways to improve maneuverability. 
Time-optimal transition to 8 = 70° reveals however an apparent singular control arc. This phenomenon 
has been generally observed for optimization of maneuvers with high values of angle of attack. 
Calculations performed so far suggest that in such situations time-optimal transition to high value of a 
or 8 does not necessarily require bang-bang type of control, at least in the final part of the maneuver. 
Normal acceleration plots as well as pitch rate show that those maneuvers are within the desired range. 

A series of optimizations of integral quadratic performance criterion (3.12) has also been executed 
in different variants. The controlled variable was either a or 8, which corresponds to either = p 3 = 
p 4 = 0 or pj = p 2 = p 4 = 0. Attempts to simultaneously regulate values of two or more states were 
not successful due to the control (particularly thrust) limitation. The simulations were performed with 
initial conditions corresponding to trim flight at T M = 3,000, 10,000 and 18,000 pounds. In this case, 
unlike in time-optimal problem, only optimal stabilator/elevator and thrust direction command signals are 
calculated. The thrust magnitude is assumed to be increased (or decreased) by the pilot manually to the 
prescribed terminal trim. The control horizon is 8s — a value which, if significantly increased or 
decreased, does not affect the shape of optimal trajectories. Typical results for control of a are shown 
on Figs. 3.6-3.13 and for control of 8 on Figs. 3.14-3.17. It is seen that the control is bang-bang in the 
initial segment of the trajectory but then smooths out and does not hit the bounds. The ripples on control 
trajectories are a result of finite discretization time A. In the numerical optimization variable A was used 
to speed-up calculations. If a decrease of A caused no significant decrease in the performance criterion, 
larger A was used. Even though the truly optimal control trajectories are much smoother than those 
obtained with discretized control, the resulting difference in state trajectories was almost unnoticeable. 
Since the controlled variable in all cases locks exactly on the desired value increasing the horizon in a 
very large range gives a picture of "semi infinite horizon control." 

Figures 3. 6-3. 9 demonstrate results of quadratic-optimal regulation of a with initial trim state for 
a 0 = 4.9° (15,000 ft altitude and 0.3 mach), and desired value of a ranging over 20.8°, 50°, 60° and 
70°. Rise time ranges from about Is to about 1.8s with a maximum overshoot less than 10%, after 
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which the control holds a exactly at its desired value. The maximum normal acceleration peaks at about 
2.5g for the most demanding case (c/ = 70°) with pitch rate peaking at less than 70°/s in about Is. 
Figures 3.10-3.13 demonstrate typical optimal maneuvers for initial conditions corresponding to thrusts 
of 10,000 and 18,000 lbs (with angle of attack at trim equal to 20.8° and 35.2°, respectively) with 
different values of target a. Since they start with larger thrust value, consequently the time needed to 
reach the desired output level is shorter. 

Optimization results for pitch angle regulation were quite similar. Again initial portions of 
optimal trajectories are of bang-bang type, and closely resemble the time-optimal case. Figures 3. 14-3.16 
demonstrate typical results for target 0 f = 60° and 0 f = 80°. Also an example of optimal descent from 
0 = 35.2° to 0=20.8° is shown on Fig. 3.17. 

Both the time-optimal and quadratic-criterion-optimal trajectories display certain similarities 
between optimal stabilator/elevator angle and thrust angle command signals, lending some credibility to 
current practice of scheduling the thrust angle proportionally to elevator angle deflection. This is 
particularly true for low angle of attack values. However, for a l> 50° these similarities begin to 
disappear. Careful examination, e.g., of Figs. 3.8 or 3.11 reveals that for truly high performance 
regulation of a at high values both control signals have to be determined independently. The same was 
found true for regulation of 0. Also in the time-optimal case Fig. 3.5 shows that both angles should be 
controlled separately. 

3.4. Practical Considerations and Conclusions 

Performed numerical optimizations provide a useful set of benchmark tests for rapid changes of 
angle of attack of pitch angle of the plane in question. Several characteristic features of optimal 
trajectories are revealed: bang-bang type command signals for actuators are required most of the time, 
with possible singular arcs for large angles of attack, even in the time-optimal case. This finding 
motivates a search for nonlinear control algorithms for rapid maneuvers, as linear controllers cannot 
successfully provide bang-bang signals. Also optimal trajectories suggest that separate control of elevator 
angle and thrust angle is advisable rather than scheduling of the latter after the former. 

For the calculations of this section a computer package was developed for solution of optimal 
control problems discussed in Section 3.2. The program was developed for the particular aircraft model, 
but due to its modular nature it is possible to substitute any plant model and any performance criterion 
of the type (3.12). Written in standard C language the package is essentially portable to any computer 
system with a C language compiler supporting ANSI standard. It is proposed that existing version of 
the program, still requiring the author's intervention to change the plant’s model, is developed into a 


3-27 



universal program allowing easy numerical solutions of optimal control problems for any plant, 
provided its model is supplied in form of appropriate C language functions. 
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4. NEURAL NETWORK CONTROL 
4.1. General Concept 

The great potential of feedforward artificial neural networks in identification, estimation and 
control of dynamical systems have long been recognized, and a multitude of applications have been 
proposed [4.4-4.6]. It is generally agreed that, although neural networks cannot be viewed as a universal 
panaceum for all design problems, they provide numerous advantages absent in other modeling 
techniques. These include automatic learning with no unnecessary assumptions about model structure, 
smooth behavior with good interpolation and extrapolation properties, and possibly very fast parallel 
hardware implementations due to simple modular structures of the networks. 

Most of the applications of feedforward networks are based on their universal approximating 
abilities. Generally, for any finite-dimensional continuous mapping there exists a multilayer feedforward 
network, with suitable size and neurons’ parameters, that approximates the given mapping arbitrarily 
closely. In particular, for control of dynamical systems, any finite dimensional feedback controller can 
be implemented by means of neural approximations. This may be advantageous if the form of the 
controller is conceptually known, but is not explicitly available in closed form. If only examples of the 
desired control actions are available, a neural network may be utilized to learn the unknown underlying 
principle from the examples. In particular, it has been suggested in [4.10] that this idea be used for 
approximation of optimal feedback using optimal trajectories calculated in open-loop mode. 

As remarked in the previous chapter, truly high-performance control calls for some kind of 
optimality. Unfortunately, optimal control theory offers mostly tools for calculation of open-loop 
controls. While closed loop solutions exist for special cases of linear systems, the design of optimal 
feedback controllers in the general nonlinear case is most often untractable. In fact, even the very 
existence of optimal feedback mapping is far from being obvious for a nonlinear plant [4.5,4. 8]. The 
approach used here is to train a neural network on a set of optimal trajectories derived numerically from 
the model of the system, under the tacit assumption that the plant is regular enough for the closed-loop 
synthesis to be indeed possible. The method from [4. 10] is extended here in the sense that a number of 
op timal trajectories is calculated not only for different initial conditions, but also for different target 
values of the controlled output variable. The numerically obtained optimal trajectories contain 
information on how the control signal depends on state variables and on the desired final state. An 
artificial neural network is trained to extract this information from the optimal trajectories and then used 
in a feedback scheme to generate a sub-optimal policy. The whole synthesis is done off-line and, 
therefore, problems with slow learning of large neural networks are not so critical here. 
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4.2. Methodology 

The neural network controller was trained to approximate quadratic-optimal regulation of angle- 
of-attack of the fourth order model of longitudinal dynamics of the modified F18 (HARV) and with all 
actuators’ dy nami cs and control constraints. The following subsections give detailed descriptions of the 
algorithms used. 

4.2.1. Feedforward Neural Networks 

The basic element of a multilayered perceptron, or a feedforword neural network, is a multi- 
input, single output static processing element, called an artificial neuron. It performs a nonlinear 
transformation: 

x -» y(x) = a (£ j wj xj + 0) ( 41 ) 

where wj are the neuron’s input weights, 8 is the activation threshold, and a is the activation function. 
By extending the input vector x with a fictitious element always equal to 1, it is possible to treat the 
threshold 8 as a weight associated with a constant input, therefore both wj and 8 are often denoted as the 
neuron’s weights. The activation function is usually a sigmoid, i.e. is smooth, monotone increasing and 
bounded, with a typical example being <r(x) = 1/(1 +e' x ). 

A network consists of neurons arranged in two or more "layers," with each neuron of the next 
layer having outputs of all neurons of the previous layer as its inputs. The neurons in each layer have 
the same numb er of inputs and outputs. Customarily all the layers besides the output one are called the 
hidden layers. The numbers of the neurons in each layer may be different and determine the properties 
of the network. The general structure of a network with two hidden layers is depicted on Fig. 4.1. 



Figure 4.1. Structure of a network with two hidden layers. 
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A feedforward multilayer network with N inputs, M outputs, K layers and with 1; neurons in the 1-th 
hidden layer can be thought of as a nonlinear function f:R N -» R M , with its i-th output equal 


„ _ JK) A K) + ^ (K) (K-l) / fl 0) . v' , 

yi - ff i 0 ij + 22j w ij a j \- e m + 22n w m,n x n- 


The properties of this mapping are determined by the weight vector 


_ _,(1) _.(1) W (K) (K) Al ) «(K) 

w- w 14 ,...,w 11N ,...,w 1>1 ,...,w M>1K _ 1 ,e 1 ,—,0 M 


which contains all neurons’ weights and thresholds. 


4.2.2. Neural Network Approximation of Nonlinear Functions 

Most of the applications of feedforward artificial neural networks are based on the approximation 
result, which may be stated as follows. Consider a compact subset of R n , U and a continuous function 
f: D -» R m . The network architecture considered here is a single hidden layer network with activation 
function a of the hidden layer nonconstant and bounded (which is satisfied by sigmoidal functions), and 
the activation function of the output layer linear. Then for any e > 0, there exists a (sufficiently large) 
network, such that the mapping (4.2) realized by the network approximates the function f in the uniform 


sense with error smaller than c, i.e. 


sup xeU 


|f° et (x) - 


f(x)| < e . 


This result is stated and proven in [4.2]. The result also holds with the output activation function 
monotone increasing and bounded with range (0,1) if the range of approximated function is contained in 
m-dimensional rectangle [0,l] m . Since a continuous function is always bounded over a compact set, this 
means that a continuous function can also be approximated by a network with sigmoidal outputs if only 
proper scaling is used. S imil ar result holds for approximation in ££? sense — i.e., for 1 < p < oo and 
e > 0 if f e i£P(U,R m ) then there exists a network such that 


| u l(f(x)-f oet (x))l p dx 


< € . 


Several extensions to this have been introduced. For example approximations in Sobolev spaces 
were also considered. If the function to be approximated is m times differentiable over U - feC m (U), and 
if the activation function of the hidden layer has bounded derivatives of order up to m, then for e > 0 
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it is possible to find a (sufficiently large) neural network approximating not only the function but also all 
its derivatives of order up to m with error smaller than e, i.e. 

E|«i<- („ l D “( f w 

or 

max | a | < m su PxeU | D “( f ( x ) " f net( x >) | < e if P = 00 • (4 ' ?) 

The assumption of compact set U can be replaced for p < co if the integration is done with respect to 
some finite measure /*. 

Approximation results reviewed above are, unfortunately, of purely existential nature, and give 
no constructive procedure to find a required approximating network. In particular, size of the network 
required to obtain desired accuracy is unknown. Also the way to determine appropriate values of the 
network’s parameters (the weights of the neurons) is not addressed. Selection of the size of the network 
is usually based upon experience or trial and error procedure, although systematic procedures for treating 
this problem have been also developed [4.9], Calculation of appropriate values of the network’s weights, 
often called tr aining of the network, is discussed in some detail in Subsection 4.2.4. 

The approximation capabilities of feedforward networks facilitate their application whenever an 
anal ytic-ally unknown function has to approximated. In that sense neural networks are in fact just one of 
man y possible methods of nonlinear function fitting, together with polynomials, splines, Walsh functions, 
Fourier series, etc. What sets apart neural networks from all those methods is their specific structure — 
with many simultaneously operating, identical processing elements, each of which carries only a small 
portion of the responsibility for the overall result. As a consequence, they are particularly well suited 
for parallel, fast acting, fault resistant hardware implementations. Also, the approximating properties of 
neural network approximations are quite advantageous when compared, e.g., with polynomials. Our 
previous studies suggest that quality of interpolation by neural networks is very good if the size of the 
network is judiciously chosen. The nature of the sigmoid function allows also for extrapolation much 
more sensible than in other approximation methods. 

4.2.3. Neural Network Sub-Optimal Feedback Synthesis 

Since a suitably large feedforward neural network can approximate any continuous function, in 
particular it can be used for implementation of a nonlinear control law. A typical situation is when the 
control algorithm, which is known to exist, is not given explicitly, but only in form of desired input- 
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output behavior of the plant, or as a set of examples of "correct" control actions. In an iterative process 
of "learning" the weights of the neurons are adjusted so that some performance criterion is minimized. 
A wealth of examples of such applications may be found in a survey paper [4.3], In this study a 
particular control application of neural networks is taken from [4.10]. 

Consider a finite-dimensional time-invariant dynamical system of the form 

dx/dt = f(x(t),u(f)) 


with control signal constrained by 


u min — — u max • 


(4.9) 


Let the control objective be either time-optimal transfer of the plants state from some initial condition Xq 
to a desired terminal state Xf, or minimiz ation of infinite horizon quadratic error criterion: 


J = \ 0 °°Ya Pi( x i' x if) 2dt 


(4.10) 


For the time-optimal control to be well posed Xf must be attainable from Xq in finite time by some control 
satisfying (4.9). Then the sufficient conditions for the time-optimal control to exist can be found in [4. 1], 
For the quadratic criterion (4. 10) minimization to be well posed, Xf should be an equilibrium of (4.8). 
In both cases, if the solution exists, it is time-invariant, i.e., if the starting time is shifted but the initial 
condition is the same, the optimal control signal is identical modulo time shift. Now consider a family 
of such optimal control problems with fixed desired state Xf and initial condition varying within some 
neighborhood of x f . This may correspond to the problem of time-optimal or quadratic optimal rejection 
of random state disturbance with fixed desired steady state value. Intuition based on the principle of 
optimality suggests that the value of optimal control signal should be expressible as a function of the state, 
hopefully with some smoothness properties. Unfortunately this does not have to be the case. Discussion 
of difficulties arising in this problem and further references may be found in [4.4] or in survey paper 
[4.8]. In this study it is assumed that indeed in the region of interest the optimal control signal for fixed 
desired state may be expressed as 

u°P‘(t) =g(x(t)) (4U) 

which is unique up to a set of measure zero (e.g., a switching hypersurface). In extension to the work 
of [4. 10] we furthermore consider a family of such feedback synthesis problems, parametrized by desired 
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final state Xf. We introduce the command or target signal Xf as additional input to the controller and look 
for a feedback law of the form: 

u°P l (t) = g(x(t),x f ). ( 412 ) 

The approach is to teach a neural network to approximate this finite-dimensional static mapping from off- 
line generated optimal trajectories corresponding to different initial conditions and different target points. 
The previously mentioned approximation theorem guarantees uniform approximation of a continuous 
function, whereas optimal controls are often bang-bang — both in time-optimal and quadratic optimal 
cases (see previous section). However, continuous functions over a compact set are dense in i£ 2 , so it 
is possible to approximate (4.12) by a neural network arbitrarily closely. 

4.2.4. Training of Neural Network Controller 

The process of finding the proper weights for the neurons is called training. It is an iterative 
gradient-based procedure, during which the quadratic error is minimized: 

' (413) 

where the sum is calculated over the training set consisting of pairs {i v yj}, with z and y being the input 
and desired output of the network, respectively, and y 1161 is the mapping realized by the network. In case 
of quadratic regulation of angle of attack, investigated in Section 4.3, the input consists of the state of 
the aircraft, and of desired value of a: z = [a, q, 6, v, af] T , and the network output is the aircraft 
control vector y = [Sj,, t m> St] t . The values of y; and z i in this case come from sampling of desired 
input-output trajectories. It has to be noted that criterion (4.13) is a discretized version of integral 
criterion (4.5) with p = 2, and this introduces yet another level of approximation. For the resulting 
network to approximate the optimal feedback mapping closely, not only has the minimum of (4.13) be 
achieved (which in itself may be a difficult and time-consuming problem), but also the training set has 
to be constructed judiciously. An obvious conflict here is between quality of the controller and effort 
required for its parameter identification (training). 

For the minimization of (4.13) the gradient of J with respect to the network’s weights is 
calculated through a chain rule, whose particular application to neural networks is nicknamed 
"backpropagation" (see [4.3]). After the gradient dJ/dw is calculated the search direction 8w (in the 
weight space) is determined and the minimum of J with respect to the step size y in the direction 6w is 
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found, which includes calculation of J for several values of y. Finally the modified values of the weights 
are set according to: 

wnew = w oW + A w = W old + y ^ 1 Sw (414) 

The flow of information during the training process is shown on Fig. 4.2. 
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Figure 4.2. Training of neural suboptimal controller. 

In this study a version of conjugate gradient algorithm from [4.7] is used for determination of 
the search direction. The line search is performed using a combined two-point gradient-based and three- 
point non-gradient parabolic fit. 

4.3. Simulations 

A neural network approximation of the optimal control was developed for the aircraft model 
discussed in Section 3.2.1. The control problem considered was minimization of quadratic performance 
index (3.12-3.13) with p 2 = p 3 = p 4 = 0, = 1 — i.e, with the angle of attack being the regulated 

output variable. A family of solutions to this problem is discussed in Section 3.3. The training data for 
the neural network was obtained by sampling 15 numerically found trajectories — with initial conditions 
trimmed with angle of attack approximately equal 5°, 20° and 35°, and with target values of a equal 5°, 
20°, 35°, 50°, 60° and 70°. The trajectories were sampled every 0.1s, which with 8s horizon resulted 
in 1200 training points necessary for a single evaluation of performance criterion. The network’s inputs 
are taken to be four states a, q, 6, v and the desired value and the outputs are Ug h , Ujj^, u^>. 
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Figure 4.3. Neural controller configuration 

A number of single hidden layer networks was trained on this training data in configuration from 
Fig. 4.2. The networks used had one hidden layer and the neuron’s activation function used was standard 
sigmoid. Then the resulting neural network controllers were tested in feedback configuration depicted 
on Figs. 4.3. A typical performance of a simulated neural network controller is shown on Figs. 4. 4-4. 7. 
Here the same initial conditions and target values were used for two neural controllers — with 30 (A) and 
15 (B) hidden neurons. The initial transition to the vicinity of required value is of very good quality. 
Then however a "steady state" error occurs. This is may be due to the fact that the network actually 
wants to perform the plant’s dynamics inversion which obviously is not 100% correct. There is no 
explicit dynamic error feedback and as a result a small error will always occur. Also the network was 
trained as a static mapping, and the approximation error minimized during training does not correspond 
to the regulation error. A much better approach would be to train the networks in the feedback loop with 
the aircraft model, accounting for different sensitivities of the output with respect to control actions in 
different state-space regions. This is proposed for our future research. 

An interesting test of the capabilities of neural controller is depicted on Fig. 4.8. Here the 
desired value of the output is set to 80° — a value that is outside the range covered by the training data. 
A quite satisfactory output trajectory adds credence to the claims about good extrapolation properties of 
neural networks. 

4.4. Conclusions and Practical Considerations 

The neural network based controller presented in the previous subsection demonstrates how the 
results of optimal control theory may be applied to the synthesis of feedback controllers. Results are 
promising, but at least a few shortcomings are apparent. First of all the simulations reveal problems with 
terminal state accuracy. These result from inaccurate inversion of the plant dynamics obtained through 
the tr ainin g process. This indicates a need for incorporating die idea of dynamic error feedback into the 
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neural controller design. One possibility is to allow a second parallelly acting regulator, possibly a linear 
one, activated in the neighborhood of the desired value. Such additional controller would provide 
corrective action to eliminate the error. Another possibility is to establish an outer loop controller 
regulating the a”* 1 value fed to the network in order to get true desired value as output. This possible 
controller configuration is depicted on Fig. 4.9. In both cases the established linear control techniques 
may be used for design of the local controller, since it is not intended to work with large variations of 
a. 


Modified plant 



Figure 4.9. Proposed configuration with an additional controller 

The more general problem involves robustness with respect to model inaccuracy. The optimal 
trajectories, used as source of training data for the network, were generated using an exact model of the 
plant. Discrepancies between the plant and the model are sure to cause serious problems in terms of the 
controller performance. For small modelling errors the previously discussed concept of additional linear 
feedback might provide a satisfactory solution. For larger uncertainties the approaches introduced in 
Sections 5 and 6 and in [4.11] may be more useful. The latter suggests a number of controllers designed 
for different no minal models of the plant and with a hierarchical classifying network interpolating between 
the nominal controllers. Work on applying that method to control of an aircraft is currently under way 
and proposed for continuation. 

The neural networks used here were merely simulations implemented on a serial computer. This 
necessarily is an impediment for the processing speed, particularly if the learning phase is concerned. 
However recent advent of hardware implementations of neural networks makes practical application much 
more realistic. Faster future computing power may make on-line training with the actual airplane 
maneuver more feasible and thus alleviate some of the problems associated with simulation training. 
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5. NONLINEAR ADAPTIVE CONTROL 
5.1. General Concept 

This section discusses results of studies of a control approach based on input-output model of the 
control plant. The majority of practically used control algorithms utilize, either explicitly or implicitly, 
some model of the plant dynamics. One of the possible forms of such a model is an input-output model, 
relating applied control action to resulting observed behavior of the system. The main reason to use 
input-output, rather than state space models, is that they employ only measured quantities subsequently 
used by the controller, and therefore are more natural in control system setting. Also, very often the 
model of the plant is not given prior to the controller synthesis, and has to be identified, either off-line 
or on-line, using the available input-output data. In such a case the input-output modelling approach is 
more effective, since it has simpler model structure and results in fewer parameters to be identified. 
Furthermore, for a given plant and given input and output signals, the state space model is not unique, 
while the input-output model is — although, of course, both may be practically realized by various 
approximate modelling techniques. Obviously, the input-output modelling approach also has its 
disadvantages. The main one is that it is basically a black-box-type technique, in which the phenomena 
"inside" the plant are of no interest, as long as its response to the input is modelled correctly. If the 
dynamics of the plant is easily available from physical considerations die state space model usually can 
be constructed with no difficulties and its parameters have well understood interpretations. On the other 
hand parameters, of input-output models usually have no immediate physical interpretation. 

There are several standard input-output modelling techniques for nonlinear systems in both 
discrete and continuous time settings. They include Volterra series, Wiener series, nonlinear time series, 
neural networks, etc. — a detailed survey may be found in [5.4], In this work the time series approach 
is used. This technique is a natural extension of discrete time modelling of linear systems, known in the 
stochastic setting as auto-regressive moving average (ARMA) models. Therefore, an often used acronym 
is NARMA — for nonlinear ARMA [5.2]. The nonlinear time series expresses future values of outputs 
as a nonlinear function of a finite number of past values of output and of control. For the purpose of 
system identification this unknown nonlinear function is usually decomposed into a sum of nonlinear 
functions with parameters to be identified appearing linearly. This allows for easy application of 
parameter identification techniques from linear systems theory, although their convergence in an on-line 
identification setting in a feedback loop is a for more complicated question than in the linear case. If the 
time-series model is to be used for calculation of control action, it is also desirable that it should be easily 
solved for current value of control. Therefore, most often control appears in the model either linearly 
(or more correctly in an affine fashion) or as some strictly monotone function. 
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In the aircraft problem, the physical model of the dynamic^ is well known and is easily 
expressible in state space form. Nevertheless, there are significant reasons to look at input-output black- 
box-type modelling as an alternative approach. The main problem arises from the aerodynamic stability 
derivatives. They are complex nonlinear functions of angle of attack, Mach number and altitude. If these 
relations are entered into the state space model it appears so complicated that its usefulness for on-line 
control generation becomes quite doubtful. Furthermore, the exact form of the dependencies of stability 
derivatives on state variables is not known — even for fixed angle of attack, Mach number and altitude 
their measurement accuracy may be as low as 50% in extreme cases. Therefore, an overly complicated 
nonlinear state space model would not seem to be very useful for on-line control calculations. A much 
better choice is a simple input-output model capturing the essential dynamics of the aircraft. In this work 
the time series model used includes polynomial nonlinearities in angle of attack. 

The control algorithm used here is an adaptive, or self-tuning, non-linear model reference 
technique. The time-series input-output model is identified on-line using a recursive least squares (RLS) 
method, and one step ahead error between the predicted and reference output is minimized. 

The following sections describe in some detail identification and control algorithms used, provide 
simulations and discuss conclusions for future research and practical implementations. 

5.2. Methodology 

5.2.1. Nonlinear Time Series Model 

The nonlinear time series model considered here has general form: 

y mod (k + l) = f(y(k), y(k-l), ..., y(k-n), u(k), u(k-l), ...,u(k-m) (51) 

Where f is a nonlinear function. To facilitate easy identification of the model the form assumed in most 
practical applications is: 

y mod (k+p) = P T <Kk) (5 ' 2) 

where P is the parameter matrix, and <j>( k) = <£(y(k),y(k-l),...,y(k-n),u(k-l),u(k-2),...,u(k-r)) is a 
nonlinear function of current and past values of input and output. Linear dependence of the future output 
on unknown parameters allows one to use one of the linear regression methods for identification [5.1]. 

The plant studied here is a fourth-order nonlinear model of the longitudinal dynamics of the 
aircraft. The same model was used in studies described in Sections 3 and 4, and is discussed in detail 
in Appendix B. The difference here is that only one control signal is used — the stabilator/elevator angle, 
while thrust is assumed to point parallel to body axis and to be of constant magnitude all the time. Thus, 
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the simulated plant has four states — a, q, 0 and v, and one control 6 h . For the purpose of the control 
system synthesis, however, measurements of only two outputs are assumed — a and q. Thus in the 
identification model we have 


y = [«,q] 

P = [PcrPq] 


(5.3) 


The choice of elements of the regressors’ vector <j> is motivated by the fact that nonlinearities in the short 
period dynamics are associated with angle of attack. Also it is recognized that due to the highly nonlinear 
nature of the aircraft dynamics it is probably impossible to fit a "black-box" model describing the plant’s 
dynamics accurately in the whole range of flight condition. Instead it is more practical to fit a simple 
local approximate model, and allow its rapid adaptation according to the change of operating conditions. 
Therefore, the regressor form is chosen as: 


<£ k = [a, a 2 , a 3 , q, qa, qa 2 , qa 3 , u, ua, ua 2 , ua 3 , l] (k) . 


(5.4) 


This results in a second-order model corresponding to a currently dominant short-period behavior, but 
including nonlinearities in a, particularly significant at large values of angle of attack. 


5.2.2. Identification Algorithm 

For the on-line identification of the unknown parameters of the model (5.2)-(5.4), a recursive 
least squares (RLS) algorithm was implemented. The basic update formula for parameter estimates at 
moment k is 


p(k) = (Q(k-2)<Kk-l)) 


/ <X(k-l) + 


<Kk-l) T Q(k-2)*(k 


-d) 


(5.5) 


with covariance matrix Q being updated according to: 


Q(k-1) = (1/X(k-1)) [Q(k-2)-(Q(k-2)^(k-l)^(k-l) T Q(k-2))/ 

(\(k-l) + tf>(k-l ) T Q(k-2) <Kk-l)l) 


(5.6) 


Since two parameter vectors p a and p q are estimated, two covariance matrices Q a and Q q are also 
calculated. Due to the approximate nature of model (5.2H5.4) the currently fitted parameters are 
expected to change together with operating conditions. To allow for such changes the forgetting factor 
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X is introduced. As a precaution against unlimited growth of the covariance matrix Q at the steady state, 
which is a consequence of X < 1 when the input is not persistently exciting, a policy with variable 
forgetting factor is adopted. The formula for X(k), taken from [5.3], is: 


X(k) = 1 - e (e(k) 2 / e(k) 2 ) 


(5.7) 


where e(k) is current prediction error of the model: 


e(k-l) = y(k-l) - y mod (k-l) = y(k-l) - p(k-l) A <#>(k-l) 


(5.8) 


and e(k) is an average of e(k) over a few last samples. For two identification processes going on for p a 
and p q two values of forgetting factor X a and X q are determined with y in (5.8) denoting either a or q. 
As an additional precaution against uncontrolled growth of matrix Q, its trace is monitored and Q is reset 
to a diagonal matrix whenever the threshold value for the trace is exceeded. 


523 . Adaptive Control Algorithm 

The control methodology is based on model algorithmic control (MAC) [5.5]. The basic formula 
of the algorithm is: 


ref (k+l) = y mod (k+l) + (y(k) - y mod (k)) 


(5.9) 


where y ref is required (reference) output, generated using the command signal (in the aircraft case from 
the pilot), y is the actual output and y mod is the output value predicted by the model: 

(5.10) 


y moa (k+l) = p T <Kk) - 

Since regressors’ vector <£( k) depends on control value u(k), the equation (5.9) can be solved for u(k) 
uni ting in required output in moment k-t- 1. In the case of the aircraft control problem, the output is 
assumed to be the angle of attack a, and the model is (5.2)-(5.4). The value of control necessary to get 
the required output value is calculated as: 


u(k) = (a r -Pi a a-p 2a a 2 -P3a“ 3 -P4«q-P5 a a -P6a ( l“ 2 -P7 a < I 0!3 -Pl2j/ 


(5.11) 


+ P9a a+ Pl0a“ + Pll 


*“ 3 ) 
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with 


a T = a ref (k+l) - (ar(k) -a mod (k)) . 


(5.12) 


If the computed control exceeds the bounds for 8 h its set equal to the maximal (or minimal) allowed 
value. If the processing time of the on-board computer needed to calculate (5. 10)-(5. 12) is small enough 
in comparison with the sampling time used, then the above formula may be used in the controller. If, 
however, the calculation time cannot be neglected, only information sampled at moment k-1 can be used 
for calculation of control at time k. Therefore, the correction term in (5.9) is replaced with prediction 
error at moment k-1 and the equation becomes: 

a ref (k-l) = a mod (k+l) + (a(k-l) - a mod (k-l)) (513) 


with a mod (k+ 1) being calculated using the estimated values of angle of attack and pitch rate at moment 
k, instead of yet unavailable actual values: 


or m<xl (k+l) = pj&k) 

<£(k) = [a, a 2 , a 3 , q, Ha, Qa 2 , Qa 3 , u, ua, ua 2 , ua 3 , l] (k) . 


Estimated current values of measured outputs are calculated taking into account previous prediction 
errors: 


m = pJ<Kk-D - (of(k-l) - a mod (k-l)) 
m = P q T <^(k-l) + (qCk-l) -q mod (k-l)) . 


Also in equation (5. 1 1), which is used for the calculation of control value, a and Q are used instead of 
a and q, and (5.12) is replaced with 


a T = ot ref (k+l) - («(k-l)-a n,od (k-l)) . 


(5.16) 


The algorithm is made adaptive, or self-tuning ([5.6]), by including the identification mechanism 
described in Section 5.2.2 for parameters p a and p q . Thus the design is a certainty-equivalence one, and 
its validity depends upon the assumption that indeed the identification and calculation of control can be 
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performed separately. Here we tacitly assume legitimacy of such method for the problem being solved, 
but it has to be pointed out that this is a heuristic method justified mostly by simulation experiments. 

The generation of reference trajectory a ref is performed using the command signal and current 
output value a. The reason to include output feedback in calculation of a ref is to avoid unrealistic values 
of reference signal, leading to control values outside the region of local validity of identified model. If 
in the previous step the reference value was not attained, then the reference in the next step is less 
demanding, so that it can be attained. In this study a second-order algorithm is used for calculation of 
a**: 


y ref (k+l) = aja:(k) + a2a(k-l) + a cmd . (5-17) 

The overall flow of information in the algorithm is schematically depicted on Fig. 5. 1 


a_cmd 

► 



Figure 5.1. Flow of information in adaptive control algorithm 

Attempts of the control algorithm to solve (5.9) exactly may result in control signal changing very 
abruptly and often oscillating, while trying to track the reference trajectory. Therefore a modified one 
step ahead design was introduced aimed at smoother control trajectories. Instead of solving (5.9) the 
following cost function is minimiz ed with respect to current control value: 


j = (a mod (k+l) -or r (k+l)) + p(u(k) -u(k-l )) 2 


(5.18) 
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with a mod , a 1 as before. Minimization of (5.18) results in 


u(k) = 



■pu(k 


-1))/ (b 2 -p) 


(5.19) 


with 

a = p la or + p^a 2 + p 3a a 3 + p^Q + p 5ot & + p^a 2 + p 7o $a 3 ♦ P i2a 

b = P8a + P9a“ + Pl0a“ 2 + Plla« 3 • 

Obviously, for p = 0 (5.19) reduces to (5.16), whereas for p — oo we obtain u(k) = u(k-l) = const. 
As before, if (5. 19) results in a value outside the allowed range, the closest feasible value is used instead. 

5.3. Simulations 

The adaptive control algorithm described in Section 5.2 was simulated on the complex fourth- 
order nonlinear model of longitudinal dynamics of the HARV aircraft. Only one control signal, the 
elevator/stabilator angle, is used with both thrust magnitude and thrust direction assumed to be constant. 
The maneuvers presented here were simulated at 15,000 feet and 0.3 Mach. Initial estimates of 
parameters p a and p q are set to zero, and the simulation is initiated at trim condition corresponding to 
a = 5°. Then the adaptive controller is simulated with a cmd = 5° for 5 seconds, during which the 
controller learns current dynamics of the aircraft. This corresponds to a practical situation, when the 
control and identification is performed continuously and before any maneuver the parameter estimates will 
never be zero. The variable forgetting factor (5.7)-(5.8) is calculated with e = 0.01 and the past 
prediction errors are averaged over last 10 samples. 

Figure 5.2 displays the result of simulation of nonlinear adaptive algorithm for command signal 
first jumping to 60° and then decreasing to 30° and 5°. The output trajectory is showed together with 
reference signal, which is displayed in dotted line. For contrast. Fig. 5.3 shows simulation of the same 
algorithm, but with linear identification model used. It is observed that controller utilizing the nonlinear 
model performs much better. As seen on the figures both these results were done on a simulation model 
that did not account for saturation of change rate of the elevator angle. If this phenomenon is included 
in the plant’s dynamics both controllers sustain undesirable oscillation. This was probably caused by 
lower maneuverability with rate saturation and by inadequate modelling of saturation with the model 
(5.2)-(5.4) which includes only terms linear in control. A possible remedy is to take into consideration 
the rate limitations in calculation of control value. Also for higher maneuverability both thrust magnitude 
and thrust direction ought to be included as control signals. This will require a little more involved 
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Figure 1 . Response for Linear MAC 






calculations for minimization of (5 . 1 8), which will become a three-dimensional optimization problem with 
constraints. These problems are addressed in Sections 3, 4, 6, and 7. 
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6. BILINEAR AND LINEAR ADAPTIVE CONTROL 

6.1 Concept 

This section describes the design of an adaptive controller for a high performance highly 
maneuverable aircraft over flight regimes including very high angles of attack. The purpose for adaptive 
control is to provide a mechanism to account for unknown changes in the system dynamics that is to be 
controlled. The goal of traditional adaptive control is to use concepts from linear theory to control an 
aircraft over a highly nonlinear flight regime. Adaptive control for a small class of nonlinear and time 
varying system is investigated in [6.14,6.15,6.30,6.37]. Model reference adaptive control usually 
includes system identification. Also, a model system may generate a desired reference trajectory. Then, 
a controller uses this information to calculate a command signal such that the output of the system follows 
the reference trajectory. A block diagram of such model reference adaptive controller is shown in Figure 
6.1. Two important elements have to be developed for an effective adaptation routine. Here, a class of 
prediction models is selected first. The prediction model approximates the dynamics of the system, and 
it has parameters that can be modified by an estimator. The estimator is the second part of the 
adaptation. It estimates the values of the parameters to improve the prediction model. The simplest class 
of prediction models are those that are linear in parameters. In this project, each model of linear, 
bilinear, and nonlinear prediction models were checked for control performance. The most common 
estimation algorithm for such models is the recursive least squares algorithm to choose parameters to 
minimiz e the mean-square difference between the prediction model and actual system. The purpose of 
making the algorithm recursive is to allow for on-line identification of parameters. 

The input reference model [6.8] is an intermediate step that allows the system to follow the 
command signal while meeting a variety of design criteria (for instance, rise time, overshoot, settling 
time, etc.). The control is calculated such that the system follows the reference trajectory, and such that 
the control signal remains within its constraints. Each block of the adaptive controller is described in the 
sections that follow. 
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Figure 6.1 Block Diagram of Adaptive Control. 








6.2 System Dynamics 

The supermaneuverable aircraft dynamics described in this section is based on a modified version 
(HARV) of the F-l 8 aircraft. The controllers consist of the stabilator and thrust vectoring. The stabilator 
angle, an aerodynamic control input to the aircraft dynamics, is useful at normal flight conditions. The 
thrust vectoring is useful at high angle of attack, low dynamic pressure operating conditions, where the 
aerodynamic control effectiveness is inadequate. The aerodynamic coefficients are considered to be 
functions of the stabilator deflection as well as the angle of attack, Mach number, altitude, and pitch rate. 
The effects of leading edge flap, trailing edge flap, speed brake, landing gear, etc, are not considered. 
It is assumed that all simulations were to occur at a nominal constant altitude of 15000 feet. The 
longitudinal equations of motion are given in Appendix B. 

The input dynamics were described by three states — thrust magnitude (T), thrust 
vectoring angle (5 V ), and stabilator angle (6^). The stabilator and the thrust vectoring dynamics include 
a velocity limiter of 40° per second for the stabilator angle, and 80° per second for the thrust vectoring 
angle [6.2]. The differential equations are stated below. 

Stabilator Angle Dynamics: 


-40 







40 


(5h cmd -«h) > 


4 

3 


The magnitude of the stabilator angle is limited according to the following equation. 


-24.0° < 6h < 10.5° 


( 6 . 1 ) 


( 6 . 2 ) 
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Thrust Vectoring Dynamics: 


S v = 


-80 


-8 


30 ( 5 v cmd -A) 

80 


3 

— < /S V ~dy\ £ _ 
3 \ v «nd v / 3 

> I 


(6.3) 


The magnitude of the thrust vectoring angle is limited according to the following equation. 


-20° < 5 V < 20° 


(6.4) 


Thrust Magnitude Dynamics: 


< 65 > 

where T cmd is command signal of the magnitude of thrust. The magnitude of thrust is limited according 
to the following equation. 

0 ^ T <; 18000 lbs . ( 6 - 6 ) 


6.3 Prediction Model 

Here we are interested in simply approximating the real system in order to synthesize an effective 
adaptive control in real time. Modeling is important since the choice of model is often the first step 
toward the prediction or control of a process. An appropriately chosen model structure can greatly 
simplify the parameter estimation procedure and facilitate the design of prediction and control algorithms 
for the process. 
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6.3.1 Class of Models 


System models can be developed by two distinct methods. Analytical modeling consists of a 
systematic application of basic physical laws to system components and the interconnection of these 
components. Experimental modeling, or modeling by synthesis, is the selection of mathematical 
relationships which seem to fit observed input-output data. Experimental modeling is emphasized in this 
section. 

Experimental models can be described by state-space, difference-operators, autoregressive moving 
average models [6.15], etc. In this section, we discuss the difference-operator representation among 
several approaches. Generally, state-space models can be generated as a set of first-order difference 
equations. An alternative description is to use a high-order difference equation with an appropriate 
difference-operator representation. In order to describe these models in a succinct manner, we introduce 
the forward and backward shift operator q and q" 1 . If y(t) denotes the value of the sequence {y(t)} at 
time t, where t E {0,1,...}, then qy(t) denotes the value of the sequence at time (t+1), and q 1 y(t) 
denotes the value of the sequence at time (t-1). That is, 

qy(t) = y(t+l) for t ^ 0 


q _1 (t) = y(t-l) for t ^ 1 ; q -1 y(0) = o 

and consequently, 

<1 1 y(t) = y(t+i) for t 5: 0 


( 6 . 8 ) 


(6.9) 


q~ l y(t) = y(t-i) for t > i ; q -1 y(0) = 0 for 0 < t < i (•) 

The first approach is to simply assume that the model can be adequately described by a linear time- 
varying system. 
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A linear time varying system can be described by the equation 


A(q"\t)y(t) = B(q _1 ,t)u(t) (6 ‘ U) 

where A(q _1 ,t) and B(q' 1 ,t) are time varying polynomials of q" 1 . A(q _1 ,t), without loss of generality, is 
assumed to be monic. Thus, A(q"*,t) could be described by the equation below. 

A(q -1 ,t) = 1 + aj(t) q _1 + a 2 (t)q" 2 + a 3 (t) q " 3 + ... + a n (t)q _n (6.12) 

This leads to a simple prediction model with the following form: 

?(t) = $(t) T 0(t) ( 613 ) 

$(t) T = [y(t-l), y(t-2), ..., y(t-n), u(t-l), u(t-2), ..., u(t-m)] 

m T = [-a x (t), -a 2 (t), ..., a n (t), b^t), b^t), b ra (t)] ( 6 - 15 ) 

The second approach is to simply assume that the model can be adequately described by a bilinear time- 
varying system. 

A bilinear time-varying system can be described by the equation below [6-7,6-29], 

rriy m z my m z 

y(t) = £ a i y(t-i) + E bj u(t-i) + X) c jj y(t-i) u(t-j) (616) 

i=l i=l i=l j=l 

where niy, m^, are the orders of the output and input. 

This leads to a simple prediction model with the following form: 

?(t) = $(t) T 0(t) ( 61? ) 

$(t) T = [y(t-l), y(t— 2), ..., y(t-m y ), u(t-l), u(t-2), ..., u(t-m z ), lg) 

y(t-l)u(t-l), y(t— 2) u(t-l), ..., y(t-m y )u(t-m z )] 
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The third approach is to simply assume that the model can be adequately described by a nonlinear 
time-varying system. 

Input-output descriptions that expand the current output in terms of past inputs and outputs 
provide models that represent a broad class of nonlinear systems. 

A nonlinear time varying system can be described by the equation below [6.6], 

y(t) = F G (y(t-l), y(t-2), y(t-n y ), u(t-l), u(t-2), ..., u(t-n z )) ( 6 - 19 > 

where F G (-) is some nonlinear function, is about as far as one can go in terms of specifying a general 
finite nonlinear system. Model (6.19) is referred to as a NARMAX (nonlinear ARMAX) model. 


63.2 Aircraft Prediction Model 

As can be seen from Section 6.3.1, several different approaches exist to formulate a prediction 
model for a nonlinear system. 

First, a linear aircraft prediction model for aircraft is considered as follows: 

&(t+l) = $(t) T 0(t-l) (620) 


*(t) T = [o'(t-l), q(t-l), or(t-2), q(t-2), a(t-3), q(t-3), ^(t), 
6 v (t), ^(t-l), 5 v (t-l), 6 h (t-2), 3 v (t-2)] 

Second, a bilinear aircraft prediction model for aircraft is given by 


a = #(t) T 0(t-l) 


( 6 . 21 ) 


( 6 . 22 ) 


$(t) T = ja(t-l), q(t-l), a(t-2), q(t-2), a(t-3), q(t-3), ^(t), 

5 v (t), 6 h (t-l), 5 v (t-l), 6h(t-2), 6 v (t-2), 

a(t-l) ^(t), a(t-2) 5 h (t-l), a(t-3) 5 h (t-2), a(t-l) S v (t)] 
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Third, a nonlinear aircraft prediction model is considered as follows: 


a = $(t) T 0(t-l) 

$(t) T = [oi(t-2), q(t— 2), a(t-3), q(t-3), «(t-4), q(t-4), 5 h (t-l), 

6 v (t— 1), 5j,(t-2), 5 v (t-2), 5j,(t-3),5 v (t-3), a(t-2)6 h (t-l), 

a(t-3) 6j, (t-2), a(t-2) q(t-2), a(t-2) q(t-2), 6 v (t-l) a(t-2), (6.25) 

5 v (t-2) a(t-3), c^(t-2), a^(t-3), a 3 ( t-2), a 2 (t-2) q(t-2), 
a 2 (t-2)5 h (t-l), V 1 ' 1 ) ^(t" 2 )] 


633 Parameter Estimation 

Overview of the Recursive Least Squares Algorithm 

The recursive-least-squares (RLS) algorithm is the most popular on-line parameter estimation 
algor ithm which is the minimization of 

J N - ±£ X N -'[y(t) -» T *(t)F < 6 - 26 > 

N t=l 

The problem is to obtain model parameter estimates which,in a least squares sense, minimi ze the 
difference between the actual output, y(t), and its value predicted by the model. The vector contains past 
input and output values and its dimension depends on the order of the model to be estimated. This leads 
to the recursive least squares algorithm with a variable forgetting factor. The parameter vector update 
law is given by 

0(t) = $(t-l) + K(t)[y(t)-d(t-l) T <Kt)] 

and the gain update by 

K( t) = 

X + <Kt) T P(t-l)<K0 


(6.27) 


(6.28) 
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Covariance matrix update: 


P(t) = 


J. 

X 


P(t-1) 


P(t-1) 0(t) 0(t) T P(t-1) 
x - <>(t) T P(t-i) <Kt) 


(6.29) 


The basic RLS algorithm with X = 1 has several important properties. First the least squares 
algorithm has a fast convergence rate (exponentially fast for a linear time invariant system with proper 
excitation). Also, the stability of the RLS algorithm combined with direct and indirect adaptive control 
is well understood and many proofs have been published in this area [6.5,6.15,6.20], 

The main disadvantage with basic RLS is that the covariance matrix gradually decays to a small 
value and therefore the algorithm does not retain its adaptivity to adequately track time-varying systems. 
The covariance matrix in the RLS algorithm tends towards zero which causes the adaptation to turn off. 
This is undesirable in the case where the parameters are time varying. Several modifications have been 
made to the RLS algorithm to correct this problem. A variety of modifications are proposed in the 
literature to keep the algorithm awake. The modifications in general are of two different types. The first 
idea is the inclusion of the forgetting factor. The second type of modifications that have been proposed 
is to mani pulate the covariance matrix directly. Modification of the covariance matrix approach is 
considered in the following section. 


Constant Covariance 

An approach proposed by Goodwin in [6. 15], is to maintain a constant covariance by the addition 
of a properly scaled identity matrix. This leads to the following algorithm. 


P'(t) = [i - K(t) <£(t) T ] — 

A(tJ 


Let r = trace (P'(t)), and C 0 , Cj denote two positive constants such that Cj > C 0 . 


(6.30) 
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IF T > C 0 

Ci ~T 

P(t) = P'(t) + — I 

n 

EF r < C 0 


(6.31) 



c i ~ c o j 


(6.32) 


t m 

The algorithm ensures a constant trace of Cj, and the following bounds are placed of the eigenvalues of 


P(t). 


c i ~ Cq 

m 


< X[P(t)] < C x 


(6.33) 


Covariance Regularization 

The basic idea of updating is a combination of a covariance resetting feature and a guarantee of 
lower and upper bounds on the covariance matrix. This algorithm replaces equation (6.29) by: 


P' 


(t) 


4 - 


K(t) <Kt) 


r] P(t-l) 
J X(t) 


(6.34) 


P(t) = 



P'(t) + C 0 I 


(6.35) 


where C 0 , Cj denote two strictly positive constants such that Cj > C 0 . 

This modification maintains the following bound on the eigenvalues of the covariance. 


C 0 < X[P(t)] < Ci (6.36) 

Its performance was reasonable, but the best results were obtained by combining the matrix regularization 
with the constant covariance. This resulted in the following algorithm. 

Let r = trace (P'(t)), C 0 , Cj denote two positive constants such that Cj > Cq, and 0 < C 2 < 1 
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IF r > C 0 

p(t) = c, p'(t) + Cl — I ( 6 - 37 ) 

n 

IF r < C 0 

P(t) = rP'(t) + Cl ~ C °I ( 6 - 38 > 

r m 

One way to interpret this algorithm is that it is a combination of the constant covariance and the 
covariance resetting. 

6.4 Control Calculation 

The controller was designed to perform or meet several goals. First and most importantly the 
control values are calculated such that the angle of attack of the aircraft follows the reference model. The 
control values were also calculated such that the thrust vectoring would return to zero if it were no longer 
needed, and a certain amount of smoothness was desired for the control signals. The following cost 
function was minimized at each step. 


1 ■ 

z z L J (6.39) 

6.4.1 One Step Ahead Predictor Controller 

Let the prediction model in equation (6.13) be described by, 

a(t+l) = a6h ft) + bS v ft) + *(t) T 0(t-l) (6-40) 

“cmd v cmd 
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S v c ^ *«] 


(6.41) 


0(t-l) T = [a(t— 1) b(t-l) 0(t-l) T ] 
Taking the derivative of J with respect to the control gives 


(6.42) 


dJ 


dS h C md 


dJ 


dS 


v cmd 

+ Pa 

This control yields the following: 


- - o(t+l)] (-a) . P2[«h c J0 -\Jt-D] 

- PlMtH-l)-&(fl)](-b) ♦ P 3 [«v c Jt)-S Vcmd (t-l)] 


(6.43) 


v cmd 


®] 




P\ a +p 2 Piab 

-1 

Pl ai» +p 2 6 hcmd (t-l) 

V (t) . 


P\ ab p 1 b 2 +p 3 +p 4 


Pi b i; + p 3 1) 


where 


(6.44) 


v = a ref -<Kt) T 0(t-D (645) 

To include the velocity and magnitude limits in the control calculation two extra conditions are added. 
The first condition requires that 6 Vcmd (t) be recalculated if \. md (t) has reached the magnitude limit. The 
second condition requires that 5 Vcmd (t) be recalculated if ^h cmd (t) is a value requiring 80° per second. 
S„ ,(t) is recalculated as follows: 

v cmd v 7 


5 v c J» - 




(6.46) 
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After the control values have been calculated they are limited by 40° per second for ^^(t) and by 80 
per second for 5 Vcmd (t). 


6.5 COMPUTATION 

6.5.1 Models 

In design of real systems, there exist some restrictions to be considered because of physical 
properties. For example, the input dynamics are described by three states-thrust magnitude, thrust 
vectoring angle, and stabilator angle. All of these have the limitations noted in Section 6.2.2. 

The linear prediction model is rewritten as follows: 

&(t+l) = b 0 (t-l)5 hcmd (t) + bjCt-1) 5 Vcind (t) + d(l-l) 

$(t) T = [ctft-l), q(t — 1 ) , a(t-2), q(t-2), a(t-3), q(t-3), 

«(t-l) T - [b 0 (t-l) bid-1) «d-l) T ] 

«‘> T ■ K s v, j» *® T 

The bilinear prediction model becomes 

«(t + l) = bud-D^O) - b 12 (t-l)6h cmd (t)«(t-l) - b 13 (t-l)6 Vcmd (t) 

«■ b 14 (t-l)6 Vcmd (t)a(t-l) + i(t) T 0(t-l) 


(6.47) 

(6.48) 

(6.49) 

(6.50) 
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(6.52) 


$(t) T = [a(t-l), q(t-l), a(t-2), q(t-2), a(t-3), q(t-3), 

5v c«n/ t_1 )’ ^cmd^ -2 ^ ^cm/*" 2 ^ 

0(t-l) T = [b u (t-l) b 12 (t-l) b 13 (t-l) b 14 (t-l) 0(t-l) T ] 


<Kt) T = 


du (t) Su 

‘ 1 cma ’ h cmd 


(t)a(t-l) S^Jt) 5,^(0 o(t-l) tft) 1 


The nonlinear prediction model becomes 


&(t + l) = c n (t-l)6 hemd (t) + + c^t-D^JOo^a-l) 

c 14 (t-l)Sv emd (t) + c 15 (t-l)6 Vcmd (t)a(t-l) + c 16 (t-l) 6 Vcmd (t) o^Ct-l) *(t) T 0(t- 


i(t) T = [a(t-l), q(t-l), a(t-2), q(t-2), a(t-3), q(t-3), ^(t-l), 

5 v (t-l), 5 h (t-2), 5 v (t-2), a(t-l) a(t-2) ^(t-l), 

a(t-l) q(t-l), a(t-2) q(t-2), 8 v (t-l) a(t-2) 

5 v (t-2) a(t-3), c?(t-2), o?(t-3), a 3 (t-2), a 2 (t-2) q(t-2)] 
8(t-l) T = [c n (t-l) c 12 (t-l) c 13 (t-l) c 14 (t-l) c 15 (t-l) c 16 (t-l) 0(t-l) T ] 


= [\J» 5 h cmd « 2 (t-D 5v cmd (t) «v emd (t)«(t-D 6 Vcmd «2(t-l) i(t) T 


6.5.2 Identification and Control 

This adaptation was performed using the modified RLS described in Section 6.2. 3. 3. 
The controller is calculated to minimize the cost function in equation (6.39) by 


(6.53) 

(6.54) 

(6.55) 

(6.56) 

(6.57) 

(6.58) 
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(6.59) 




Pi a 2 +P 2 

Pi ab 

-1 

Pia 1 , + p 2 6 hcindt -l 

V (t) . 


Pi ab 

Pi b 2 +p 3 +p 4 


Pl bit + P 3 $ Vcmd (t-l) 


This leads to the following control law calculation: 

5 hcmd (t-l)(b 2 / : > 1 P2 +P2P3 + P2Pa\ + (ai/ -ab 5 Vcmd (t-l)piP 3 + a i?p 1 p 4 
s h ca Jti = — 5 — 5 5 — 

b pj p 2 + a p { (p 3 + p 4 ) + p 2 (p 3 + p 4 ) 


= 


_ 6 v cind ( t - 1 )( b2 PiP 3 + P 2 P 3 ) + ( b »?-ab6h cind ( t - 1 )PiP2 


(6.61) 


b 2 p 2 p 2 + a 2 p 1 (p 3 +p 4 ) + p 2 (p 3 + P 4 ) 


where 


V = «ref ~ 


(6.62) 


6.6 Simulation 

Several different simulations were used to evaluate the model performance. Two cases of the 
results are described below. The first maneuver corresponds to the maneuver presented by Ostroff in 
[6.32], The angle of attack is changed from 5°, to 60°, 35°, and back to 5° in 8 second intervals. In 
the second maneuver the angle of attack is changed rapidly from 5° to 85° for an extended period of 
time. 


6.6.1 Simulation Data 

In the case of the linear prediction model, the controlled maneuvers use C 0 = 600, = 1200, 

Cj = 0.98, p l = 100, P2 = 0.001, p 3 = 0.001, p 4 = 0.001. 
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In the bilinear prediction model, the first maneuver starts with C 0 = 800, C L - 1600, 
&2 = 0.98, p : = 92, p 2 = 0.0001, p 3 = 0.0099, p 4 = 0.0. 

For the second maneuver, the bilinear prediction model starts with Cq = 800, = 1600, 

C 2 = 0.98, p x = 89, p 2 = 0.07, p 3 = 0.00001, p 4 = 0.00044. 

In the nonlinear prediction model, both maneuvers start with C 0 = 1200, C x = 2400, 
C2 = 0.98, Pl = 95, p 2 = 0.0001, p 3 = 0.01, p 4 = 0.0. 

6.7 Conclusions 

The character of the response for the first maneuver in the linear prediction model, the bilinear 
prediction model, and the nonlinear prediction model is similar to the response reported by Ostroff in 
[6.32]. However, the adaptive controller provided a somewhat faster response. The angle of attack 
reached 55° in approximately 2.0 seconds while the variable gain approach in [6.31,6.32] reached 55° 
in just under 3.5 seconds. The time optimal control (with a limitation of 40° per second on the thrust 
vectoring) reached 55° in about 1.8 seconds [6.24]. In the case of maneuvers one and two, comparing 
the bilinear adaptive controller with the linear controller, we see that the response of the bilinear adaptive 
controller is slightly faster for given desired trajectories (60°, 35°, and 5°) and has a smaller oscillation 
at the terminal end. For the second maneuver, the bilinear controller held the angle of attack at 85° as 
the aircraft approached steady state. 
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7. SLIDING MODE CONTROL 

7.1 Concepts 

Sliding-mode control (SMC) appeared in the Russian literature as early as 1930s and was further 
developed independently in Russia and the U.S.A. along with variable structure control (VSC) systems 
in the early sixties (e.g., [7-30,7-49]). VSC systems is a part of nonlinear control systems where the 
structure is not fixed but is varied as part of the control process, for example, by switching control-gain 
between two values according to some law. 

The SMC typically uses a high-speed switching control law to make the state trajectory of the 
system approach a switching ("discontinuous”) surface which is called the sliding surface and to remain 
on this sliding surface once the state trajectory intersects the surface, or at least for some desired period 
of time such as in handling state constraints. 

Design of SMC consists of at least two main parts. One is the reaching mode, in which the state 
trajectory starting from anywhere on the state plane, is directed toward the sliding surface. The other 
part is the sliding-mode with the state trajectory on the sliding surface. For these two modes, we have 
to design the switching surface and control law for reachability and existence of a sliding-mode. 

In general, discontinuous control is used and it is switched between two values near the sliding 
surface to satisfy the reachability condition. In the ideal case, the switching time of a sliding-mode is 
zero, but in reality fast switching of the control input causes the state trajectory to chatter along the 
sliding surface. This imperfection in the switching mechanism not only generates the undesirable high- 
frequency components in the state trajectory but also may excite unmodeled high-frequency system 
dynamics. Consequently, it could make the system unstable, and in many cases, chattering must be 
minimized. To reduce the effect of chattering by smoothing out discontinuous control, various 
investigators (e.g., see [7-46]) introduced a boundary layer near the sliding surface. The boundary is 
defined by the set L(x) = {x | ||s(x) || ^ e, e > 0} where e is the thickness of the boundary layer and 
s is sliding surface. 

Outside the boundary layer, control input is chosen such that attractiveness to the boundary layer 
is guaranteed. Inside the boundary, discontinuous input sgn(s) is replaced by (s(x)}/e. 

Generally, there is a trade off between a fast response and overshoot in control design. However, 
SMC potentially can alleviate this trade off. If states are driven to a sliding surface and slide along it 
keeping the angle of attack at desired value as in Fig. 7-1, the state finally arrives at a desired point p 
without overshoot in an ideal sliding-mode control. 
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Fig. 7.1 Ideal sliding-mode control 

If continuous sliding-mode control (CSMC) is used and the reachable speed is adjusted within 
a boundary layer the effect of chattering can be reduced and a better result can be obtained without 
overshoot. This nonlinear control law uses a sliding surface which is described in the form of the 
equation of error between actual output and desired output. This control law is then applied to the 
nonlinear longitudinal motion of a high performance and highly nonlinear aircraft, HARV (F-18), and 
as a result, the rise time can be reduced and the output reaches its final terminal point without overshoot. 
CSMC also shows some robustness to parameter uncertainty and to disturbances with given bounds. 

7.2 Methodology 

Usually, sliding-mode control with discontinuous input which includes the sgn( • ) function induces 
chattering. To reduce the effect of chattering, a continuous control law is adopted. This continuous 
control law can be expanded to not only the decoupled input-output system but also to a multi-input multi- 
output system. Further more, it can be applied to systems with constraints of control by modification. 

The main purpose of this control is for rapid maneuvers with large changes in angle of attack and 
to keep it at the final value while the aircraft is not trimmed. Here, only two controlled outputs — angle 
of attack a and pitch rate q will be selected. During this control, the thrust magnitude will be scheduled. 
The desired trajectories of angle of attack a and pitch rate q are selected so that pitch rate converges 
when angle of attack approaches its final terminal. The control scheme is composed of a-q control for 
rapid variation of a and a-control for settling a at its terminal value. 
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The control scheme will be described in detail for a difficult but representative maneuver for 
increasing(or decreasing) angle of attack rapidly between 5° and 60°. Of course, the same control, with 
appropriate simple algorithmic adjustment, is just as successful for general classes of maneuvers. 


7.2.1 System Dynamics 

The fifth-order longitudinal motion of the HARV (F-18) described in Appendix B can be 
represented simply by the state model: 


« = f l + gll 5 e + gl2 T x + gl3 T z 

(7.1) 

q = h + g21 5 e + g22 T x + g23 T z 

(7.2) 

v = f 3 + g 3 i 5 e + g 32 T x + g 33 T z 

(7.3) 

e = q 

(7.4) 

h = V sin(y) 

(7.5) 

where 

initial value [a 0 q 0 V c 0 O hj = [4.3° 0 4.3° 500 15000], 



a : angle of attack (degree), q : pitch rate (degree/second), V : total speed (feet/sec), 

6 : pitch angle (degree), 6 h : stabilator angle (degree), y : climb angle (= 6-a), 
h : altitude (feet), T x : x-direction thrust magnitude (lbs), T z : z-direction thrust magnitude (lbs), 
and <5 V : thrust angle (= tan^OyT*)) (degree). 

F = [fj f 2 f 3 ] T < 7 - 6 ) 

and 


G = [Gj G 2 G 3 ] = 


gu gl2 gl3 
g21 g22 g23 
g31 g32 g33 


(7.7) 


are functions of angle of attack, pitch rate, speed, mach number and altitude through the corresponding 
stability derivatives. 
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Assumed limits of the stabilator angle and the thrust angle include those used in Section 3 . Here, 
however, the admissible thrust angle constrained is 80°/sec. Since this was the most recent information 

available. 

7.2.2 Design of Controller 

While the control algorithm can be computered for any given flight path, an initial a 0 = 4.3 

and terminal a = 60° is discussed as typical end conditions. 

The purpose of this CSMC is to drive the concerned states to a sliding-surface which is 
represented by errors between the concerned states and their desired states. For this, control inputs <5^ 
5 V , scheduled-thrust magnitude are used and all the states (a,q,0,V,h) are feedbacked as shown in Fig. 
7.2. The control scheme is composed of two parts. In the first stage, an a-q control scheme is used until 
the error e between the final terminal value and a is less than some error e assumed to be 5 , that 
is, a arrives at nearly 55° while a and q follow the desired a d and q<j. During this stage, stabilator 
angle, thrust angle and scheduled-thrust magnitude are used as control inputs. 

Obviously, there is some trade off between a and q. Here, the major state-variable is angle of 
atts.Hr not pitch rate. Therefore, the a-q control scheme is changed to an a-control scheme when the 
error e between and a is less than 5°. In this stage, a fast approach toward the sliding surface may 
bring out a big oscillation because of constraints of control inputs. Therefore, the reachable speed should 
be adjusted using error e = a d -a and time derivative of angle of attack, a, to make angle of attack 
approach its sliding surface slowly. During this stage, stabilator angle, thrust angle and scheduled-thrust 
magnitude are also used as control inputs. 


Scheduled 



Figure. 7.2 Control structure of CSMC 
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g-q Control 

Select an error sliding surface 

s - [(“ - a d) (q - <id) (v - v d )] - ° 


and the time derivative of it assuming that V = V d is given by 


s = (a- d d )(q - q d )(V - V d )] = [s^O] =0 


Algebraically, s can be separated as follows with S3 = 0 

s = (6f - df d ) = (f t -d d + g n 5 e + g 12 T x + g 13 T x ) 

S 2 = q-Qd = ( f 2-qd + S21 5 e + g22 T x + g23 T x) • 
The reachability condition takes the following expression: 


• T 
ss 


s(F' + Gu) = sF' + sGu = sG 


u + 


(sG) T , p . 

UsG | 2 


where 

U = [Oj, T z T x | t , 


is Euclidean vector norm. 


G = [G x G 2 G 3 ] = 


gll gl2 gl3 
g21 S22 S23 
g31 g32 g33 


and 



” /* 
Fl 


( f l - “d) 

F' = 

*2 

= 

f 2-Qd 


*3 


_f 3 -v d _ 


(7.8) 


(7.9) 


(7.10) 

(7.11) 


(7.12) 


(7.13) 
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If 


u = - sF' - K(sG) t ( 714 ) 

II sG « 2 

then ss T = -(sG)K(Sg) T < 0 where 

kj 0 0 

K = 0 k 2 0 

0 0 k 3 

is a positive definite matrix. 

But ss T can be modified for practical application with bounded control inputs such as 

ss T - [M, M 2 M 3 ] fu * gli SF' - M[u +L] P.15) 

Mi + M 2 + M 3 

where 

SG = M = [Mj M 2 M 3 ] 

II M || 2 = M 2 + M 2 + M 3 2 
and 


( 7 . 18 ) 


ss T is a function of 5 h and T z because T 2 = T 2 + T 2 



( 7 . 16 ) 

( 7 . 17 ) 
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By equation (7. 15) ~ (7. 18), the reachability equation for the system (7.1)~(7.5) can be expressed by 

- Ml(«h * L,) * Mj(T z * Lj) * Mjfr, * m 

= + M 2 (T z + L 2 ) + M 3 

= s Sl($h) + 

That is, ssj is a function of <5 h and s $2 is a function of T z or 6 V . 

The control scheme is composed of two parts to satisfy the sliding-surface reachability condition 
ss^) < 0 and < 0 separately. 

At every step where sampling time At = 0.03, stabilator angle 6^ is constrained to the interval 
t^hmin’ 5 hmax] 311(1 T z is constrained to the interval [T zmin ,T zmax ] where 6 V = tan 1 ^ /T x ). 

(1) Control for ssj (6 h ) < 0 

For ss^Sj,) = + L : ), if Sj, = -L x - kjMj is selected where k t is a positive constant 

ssi(6h) < 0 is satisfied. 

Case 1: -1^ £ <5 hmax ] and ssj > 0 or ssj < 0 regardless of bounded input 6 h 

If -Lj is outside of [Shmin> ^hmaxl one of 5 hmin 311(1 ^hmax is selected to make a and q approach 
the error sliding surfaces as fast as possible considering the sign of Mj. That is, 6 hmin and 6 hmax is 
selected as inputs for Mj > 0 and Mj < 0, respectively. 

Case 2. -Lj G [5},,^, ^hnuuJ 

We choose 6 h depending on the sign of Mj 

3 h = - L j - kjMi = Shnun 2 - " — , if > 0 

and 

«h “ ^ - k i M l = 5hma |~ — > if M i < 0 

such that the reachability condition ss T < 0 is satisfied. 

After 6j, enters into the inside of the input bound, kj is calculated by 

ki = dhmhl .lh if Mi > 0 

1 2Mj 1 
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k, = 6 hmax + L 1 {f M , < 0 

1 2M, 1 

At the next step, k, is used in case that -L, is within the input bound. If k, is so large that 6 h 
is outside of [Shmin, 6 hmax ], k, is reset to some value by the above method again to make 6j, be within 
the input constraints satisfying ss < 0. 


(2) Control of ss^Oy < 0 

At every step T z = T sin($y) is bounded because 6y is constrained to l^vmin * fnr 

-(*72} < 5 V < tI 2 . 

Case 1: -1^ £ [T zmin , T zmax ] and ss 2 > 0 or < 0 regardless of T z 

One value is selected between T zmin and T zmax to make a and q approach their sliding surfaces. 


Case 2: -a < s^ < b for T z and a,b > 0 

In this case,the reachability condition is -a < ss 2 ^ 0. Therefore, the range of T z to satisfy 
-a < s^ < 0 is changed to T zl < T z < where T zl is T z which satisfies one between -a = ss 2 and 
s^ = 0 and is T z which satisfies the remained one. 

The relation of T zl , T*, T zmin and T zmax is T zmin < T zl < T z ST* < T zmax because 
T z = T sin(5 v ) is increasing function for -(ir/2) < 8 V < *72. The control input is calculated by the same 
method as that applied in case that ss, < 0 with T z G [T zl , T^] which is given by 

T z = -L2 - k 2 M 2 = Tzl ~^ , if M 2 > 0 

T z = -L2 - k 2 M 2 = Tz2 ~- , if M 2 < 0 

and different k 2 depending on the sign of M 2 which is given by 

k 2 = Tzl — if M 2 > 0 
2 2M 2 2 


k 2 = 


T z2 +L 2 

2M 2 


if M 2 < 0 


Case 3: Case that ss 2 < 0 and -L^ G [T zndn , T zmax ] 

The above same method is used with T z G [T zrnin , T zmax ]. 
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This control scheme can be explained simply in flow chart as shown in Fig. 7.3. The control 
of ssj (6 h )<0 is depicted by Loop 1 in Fig. 7.4. 


of-Control 

The a-q control scheme is switched to a-control scheme when error e = - a between final 

terminal and a is less than 5° for keeping angle of attack near its final terminal. 

According to CSMC scheme, select the sliding surface s = (a-a d ) = 0. 

The equation (7.1) can be expressed using thrust angle 

a = fi + gu 5h + g 123 sin(a-5 v ) (7.20) 

Using equation (7.20), the reachability condition is given by 


If 


ss 


Sgll 




7^4 ( fl ‘ *4 + g 123 Sin ( a " 6 v)) 

Iglll 


(7.21) 


5h = --^2 ( f l “ “d + 5123 sin (« " «v)) " k l s Sll 

Iglll 

then ss = -k^sgu) 2 <0, where kj is a positive constant. 

The equation (7.21) can not be applied directly because the control inputs - stabilator and thrust 
angle are bounded at every step as Shmin — 5h ^ «hmaxand5 vmin £ 5 V < 5 vmax . 

Let 


Then, ss = sg 11 (6 h - 1) 


l = 


- ( f i "“d + gi23 sm { a ~ 6 v)) 

Iglll 


(7.22) 

(7.23) 


Here, the control purpose is to make the angle of attack approach its sliding surface slowly 
without big oscillation satisfying the reachability condition ss<0. For this, the reachable speed will be 
reduced slowly for some time, which means that & is decreased satisfying ss < 0. To adjust reachable 
speed, multiple boundary layers are set in the neighborhood of sliding surface and the reachable speed 
is adjusted in each boundary layer satisfying the reachability condition as follows: 


• 

| error = Of ina j - a \ > 1 ° : adjust reachable 

speed 

so that 


a = ±0.1°/At 



• 

0.2° < | error = - a < 1° : adjust reachable 

speed 

so that 


a = ±0.01 °/At 
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Figure 7.3 Flow chart of a-q control 




Figure 7.4 Flow chart of Loopl 







• 0° < | error = Of- inal - a | < 0.2° : adjust reachable speed so that 

d = ±0.0001 °/At 

where At (0.03 second) is sampling time. 

When d becomes zero approximately interval [S hmin , 6hmax] enters into [/ min , / max ] or interval 
[5h m in, 5 hmax ] intersects interval [l^, / max ]. From this time, the reachable speed of angle of attack can 
be adjusted linearly in time, a-control scheme is expressed simply by the flow chart in Fig. 7.5. 

When the aircraft maneuvers from high angle of attack to low angle of attack the same method 
is applied except that angle of attack should be decreased with maximum speed by ar-q control. For 
uncertain system, CSMC also can be applied by adjusting the reachable speed in the neighborhood of the 
sliding surface. 

7.3 Simulation 

CSMC scheme is applied to the MIMO system with constrained control inputs. In this case, 
desired trajectories have to be selected appropriately considering the characteristics of the system. 

Scheduled thrust magnitude and desired trajectories of angle of attack and pitch rate for 
maneuvers from initial a = 4.3° to a = 60° and a = 60° to a = 5° have to be selected appropriately 
considering the system dynamics. Here, it is assumed that the pilot adjusts the scheduled thrust 
magnitude for the desired maneuver as shown in Fig. 7.6. |T| = T 0 = 1467.19 lbs for 0.96 second 
and | T | is increased linearly from T 0 to 18000 lbs for 0.96 < t <2.94 sec, and it is kept at 18000 lbs 
after 2.94 second. Then thrust magnitude is kept at 18000 lbs to 7.98 second. After 7.98 second it is 
decreased linearly from 18000 lbs to 10234.584 lbs for 7.98 < t < 8.93. Constant thrust 
T = 10234.584 is scheduled for 8.93 < t < 15.84 second and for 15.81 < t < 16.65 second, it is 
decreased linearly and finally reaches T = 3000 lbs. 

The main purpose of this CSMC control is to control the angle of attack as fast as possible and 
keep it at some terminal value. After the rapid a-control phase, it is assumed that the pitch rate should 
approach zero. 

The desired trajectories of angle of attack a and pitch rate q are selected as follows: 
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a-control 



Figure 7.5 Flow chart of a-control 





x!0« 



Figure 7.6 The scheduled thrust magnitude(lbs) of maneuver from ot 0 =4.3° to 


ot=60° and from a=60° to oc=5° 



Figure 7.7 The desired angle of attack(deg) of maneuver from 0 ^= 43 ° to a=60° 
and from a=60° to a=5° 



Figure 7.8 The desired pitch rate(deg/sec) of maneuver from oc 0 =4.3° to oc=60° 


and from oc=60° to a=5' 





60 -ot c 

= 2.49 * + ’ 

0 < t < 2.49 

a d = 60 , 

2.49 < t ^ 7.98 

-25 ^ 25 x 8.01 ^ ^ 

0 'j = t + + OU , 

d 8.97-7.98 8.97-7.98 

7.98 < t < 8.97 

a d = 35 , 

8.97 < t < 15.81 

-(35 - aJ 

“ d = 16.8 - 15.81 “ 15 81 ’ 35 • 

15.81 < t < 16.80 

a d = 60 , 

16.80 < t 


where a Q is the initial angle of attack. Figure 7.7 shows the desired angle of attack. 


9d 

9d 

<ld 

<ld 

<id 

9d 

9d 

9d 

9d 


50 

0.81 


t , 


-50 


t + 


50 x 3.99 


3.96-0.81 3.96-0.81 * 

0 , 

-22.5 


8.73 -7.98 
22.5 


11.61 -8.73 
0 , 


(t - 7.98) , 
(t- 11.61) , 


-17 

16.59 - 15.81 
17 

17.31 - 16.56 
0 , 


(t- 15.81) 
(t - 17.37) 


0 < t < 0.81 

0.81 < t <£ 3.96 

3.96 < t < 7.98 

7.98 < t < 8.73 

8.73 < t <; 11.61 

11.61 < t < 15.81 

15.81 < t < 16.56 

16.56 < t < 17.37 

17.37 < t 


(7.25) 


Figure 7.8 shows the desired pitch rate. 

CSMC is applied to maneuvers from low angle of attack to high angle of attack and vice versa. 
For example, from initial a = 4.3° to a = 60° and a — 60° to a = 5°. CSMC shows a fast response 
(Figs 7.6-7.13). Angle of attack arrives at 55° approximately in 2.19 second and 59° in 2.6 second 
approximately and approaches the final terminal state without overshoot. However, this control shows 
a small overshoot for a maneuver from a = 60° to a = 35°. This overshoot can be reduced by using 
appropriate desired trajectories of angle of attack and pitch rate. 
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Figure 7.9 Response of the pitch angle(deg) for maneuver from a 0 =4.3° to a=60° 
and from <x=60° to a=5° 



Figure 7.10 Response of the total velocity(feet/sec) for maneuver from a 0 =4.3° to 
a=60° and from a=60° to a=5° 



Figure 7.11 Response of the altitude(feet) of maneuver from 0t,=4.3 o to a=60° and 


from a=60° to a=5 








7.4 Conclusion 


An effective control design methodology using continuous sliding-mode control (CSMC) of a 
highly nonlinear maneuverable high performance aircraft, the HARV (F-18), has been presented. 

Although not shown here, our simulations show that with only horizontal stabilator control, 
CSMC works very well without overshoot for maneuvers up to about a = 20°. CSMC of the MEMO 
longitudinal motion of HARV is demonstrated successfully by accurate computer simulations. Similarly, 
Ostroff [7.39,7.40] investigated the maneuver by utilizing numerous trim-state linearization studies 
accompanied by scheduled variable gain in PIF controller. CSMC shows, in particular, a fast response 
from = 4.3° to 55° in 2. 19 sec. and settling time to 60° in about 3 sec. compared to the latter case, 
which shows a rise time from to 55° in about 3.3 sec. and settling time in about 6 sec. 
Accordingly, CSMC shows that it is very useful for the terminal approach to the high angle of attack 
reducing overshooting and chattering of control even though control inputs are constrained physically. 
The response using time-optimal control, of course, is faster, that is, the slightly faster as shown in 
Section 3. Therefore, better results can be obtained if time-optimal control is used for the fast response 
in the first stage, and CSMC used for driving output to the final terminal value without overshoot and 
keeping it in the neighborhood of the sliding surface. CSMC also shows some robustness to parameter 
uncertainty and disturbance with the known bounds. This is accompanied by a boundary layer in the 
neighborhood of the sliding surface and adjusting the reachable speed within it so that less chattering of 
control results for the system with constrained inputs. 
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8. CONCLUSIONS AND DIRECTIONS 

Much more was accomplished with respect to high-performance aircraft control by this project 
than was originally proposed. However, less was accomplished on stability methodology development 
than anticipated. Although, it is apparent that stability criteria, such as those desired for bilinear systems 
with nonlinear feedback [8.1], can be applied to the aircraft dynamics with feedback, we were not able 
to find relevant practical design criteria in the time allowed. Unfortunately, the desired sufficient 
conditions are overly constraining. On the other hand, a new linear perturbation study was made for a 
second-order short-period model with a range of admissible feedback gains derived in conjunction with 
that model. This result is reported in our Phase 2 Annual Report [8.2]. 

The research presented here does show the potential role of nonlinear control for high 
performance aircraft. Indeed, the time-optimal and quadratic-optimal controls investigated make full use 
of nonlinear dynamics, as shown in Section 3. Nonlinear adaptive controls (Sections 4 and 5), neural-net 
controls (Section 6), and sliding-mode controls (Section 7), all use nonlinear feedback to effectively 
control, in many cases nearly optimally, the nonlinear longitudinal aircraft motion [8.3]. 

The neural-net-based nonlinear feedback control was trained according to simulated optimal flight 
trajectories in Section 4. Unfortunately, the lengthy computations can not be accomplished on-line with 
present technology. Still, if used in conjunction with more classical on-line feedback control, such 
intelligent controls could be effective in a practical sense. Other nonlinear adaptive schemes (Sections 
5 and 6) can be done in real time, but do depend on sufficient disturbance (naturally or in terms of 
persistent self-excitation). In some cases such excitation may not be available or not tolerated and, again, 
off-line "learning" may be used. 

The tradeoffs between on-line and off-line computations, between natural and self-excitation, 
between robustness and optimal performance all need to be investigated in more detail. Convenience to 
the pilot in using the controllers also needs to be studied. 

Perhaps of immediate interest for controller, and even aircraft and engine configuration design, 
is the optimization analysis and the developed optimization software. This methodology and 
corresponding software can provide information on the most performance sensitive design parameters and 
constraints as well as providing precise control and performance data. 
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APPENDIX B 

Equations of Motion and Aerodynamic Model 



B. EQUATIONS OF MOTION AND AERODYNAMIC MODEL 
B.l Introduction and Notation 

In the following the equations of motion of an airplane in the longitudinal mode will be derived 
from the basic six degrees of freedom equations of a rigid body. The curve fitting technique for the 
stability derivatives will also be presented. 

Coordinate System 

Figure B.l depicts the body axes coordinate system used in this work. 

* 



Figure B.l. Coordinate System (body axes) 

XYZ: Aerodynamic forces in X b , Y b , directions 
FR, FQ, FP: Aerodynamic moments in X b , Y b , Z b 

The aerodynamic lift and drag forces are defined in Fig. B.l. Also, the flight path angle ( 7 ), 
pitch angle (0), and angle of attack (a) are related, as shown in Fig. B.l. 

The side slip angle /S is the angle between V T (the velocity vector) and the X b Z b plane. Angle 
of attack, a, is the angle between the projection of V t on the X b Z b plane and X b . 
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jij 


The following relations are defined 


a = tan -1 

w 

/S = sin 1 

V 


u 


V 


V 2 = u 2 + 


w 


(B.l) 


Nomenclature 


b, b' 
c 

c. g. 

c ij 

C m 

C D 

C L 

D 

F x> F y» F z 

FP, FQ, FR 

g 

h 

^yy> ^zz 


^x£> ^yf> ^z( 

P P P 
*x» *y> *z 

M 

M x , My, M 2 

m 

P 

q 

q 

r 

S 

T T T 
*x> A y> *z 


reference length for lateral derivative (span, ft) 
reference length for stability derivatives (MAC, ft) 
center of mass (gravity) 

coefficients, functions of moments of inertia, 1^, Iyy, 

pitching moment coefficient 

drag coefficient 

lift coefficient 

drag force (lb) 

external force components Ob) 
angular acceleration components Ob-ft) 
gravitational acceleration (ft/sec 2 ) 
altitude (ft) 

principal moments of inertia (slugs -ft 2 ) 
cross product of inertia (slugs-ft 2 ) 
lift force Ob) 

position vector from c.g. to engine thrust center (ft) 
position vector from c.g. to aerodynamic center (ft) 
mach number 

external moment components Ob-ft) 

aircraft mass (slug) 

roll rate (rad/sec) 

pitch rate (rad/sec) 

dynamic pressure Ob/ft 2 ) 

yaw rate (rad/sec) 

reference area ((wing area) ft 2 ) 

engine thrust components Ob) 
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u 

V 

w 

X, Y, Z 


velocity component in x b direction (fit/sec) 
velocity component in y b direction (ft/sec) 
velocity component in z b direction (ft/sec) 
force components in x^,, y b , z b directions (lb) 


a 

& 

7 

0 , <i> 

p 

5h, 6a, 6r 
a 


angle of attack (rad or degree) 
angle of side slip (rad or degree) 
flight path angle (rad or degree) 

Euler angles (pitch, yaw, roll) (rad)* 
air density (slug/ft 3 ) 

deflection angles of stabilator, ailerons, and rudder (•) 
standard atmosphere density ratio 
throttle setting (•) 


B.2 General Equations of Motion (GDOF) 

The general, GDOF equations of motion derived from Newtons laws are given in many 
references. The form used in [B.l] is shown here for the force and moment equations and in [B.2] for 
the Euler equations. 


Force Equations 

F x = m (u - vr + wq) 
52 F y = m (v - wp + ur) 
52 F z = m (w - uq + vp) 


Moment Equations 

52 M x = p - r I,^. + (Izz - lyy) q r - pq Ixz 

52 My = qlyy + rp^-I^) + Ip 2 - t 2 ) 

52 M z = fl n - pl« + (Iyy -I„)pq + Ixzqr 


*Get order of rotations: yaw, pitch, roll 
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Euler Equations 


6 = qcos <f> - rsin<£ 

<j) = p + qsin<£tg0+ rcos</>tg0 
$ = (q sin 0 + r cos </>) sec 0 

Here the aircraft is treated as a rigid, symmetrical body (so only the 1^ cross product exits); also (j mean 
d/(dt) the time derivative. 

The GDOF equations can be written in terms of u, v, w and Dyr [B.3]. 


Force Equations (body axes) 


u=rv-qw-gsin0 + X + — 

m 


v 


= pw - ru + gcos0sin<£ + Y + 



w = qu - pv + g cos 0 cos <j> + Z + 



(B.5) 


Moment Equations (body axes) 

p - C 41 pg * C 42 qr * C 43 FR * C * FP * ^((,,1, - f y ,T x ) + (< y , T, - t zl T y ) 

A XX 


*zz 


q = C 51 pr + C 52 (r 2 -p 2 ) + FQ + -L (f z/ T x - t xl T z ) 

lyy 


(B.6) 


i = C 61 pq - C 62 qr + C 63 FP + C * FR + (t yl T z - l x , T y ) + (t x£ T y - f yt T x ) 

A xx x zz 


The three Euler equations (B.4) remain unchanged. 
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In these equations 

X = [-D cos a + L sin a]/m (normal aerodynamic acceleration) 

Z = [-Dsina - Lcosa]/m (axial aerodynamic acceleration) 

Y = qSCy/m (aerodynamic side acceleration) (B.7) 

and 


L — q S Cj_ (lift) D — q S Cjj 

The aerodynamic accelerations are 

(drag) 


FP = 

[qSbCf + m (f y Z - f z Y)] / I M 

(rolling) 


FQ = 

q Sc Cm + m (f z X - f x 2)]/ 

(pitching) 

(B.8) 

FR = 

q Sb Cn + m (f x Y - f y Z)] / 

(yawing) 



Note that reference length for roll and yaw is b (wing span) and for pitch is c, mean aerodynamic chord 
(MAC). 

The inertia moments are assumed to be constant and their coefficients are defined as follows 

c* - C w = I„ • I = /(l„I„-I*) 

C„ - C*I„(I„*I„-I„)/i„I = 

C 42 * C ' [l=(l w -I4-4]/lx»Ixz 

C43 = ^ ^xz / ^xx 

/ (B.9) 

C51 ^zz-Ixxj/lyy 

C 52 = Jjjz/ 1 ; yy 

c« - c - [i„ (I„ - 1^) * 1^] / i = 

C 62 = C * *xz (*yy “ *zz " *xx) / J xx J zz 

C 63 = c *1^/1^ 
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B.3 Form of Aerodynamic Coefficients (Stability Derivatives) 

The aerodynamic coefficients, also called stability derivatives, are written in the following form: 

Drag Coefficient C n = C D (a,M,h,6h) 

Lift Coefficient C L = C L (a,M,h,Sh) + [C L (a,M,h) q + C L . (a,M,h) dl m.9) 

o 2V L <1 a J 

Pitching Moment C m = C mQ (a,M,h,5h) + [C m<i (a,M,h) q + C m ^ (a,M,h) dj 

The aerodynamic coefficients Cj^, Cl o , C mo depend on angle of attack, a, mach, M, altitude, h, and 
control surface angles, 5, (which are the controlling factor). The damping coefficients C L 0 > C La> C m a > 
C m ^ are dependent on a, M, h but not functions of 8. 


B.4 Equations of Motion for the Longitudinal Mode 

The case in which the airplane moves without side slip and rolling motion is called the 
longitudinal mode. In this the motion is restricted to a plane containing the X b Z b plane. For this case 
we may make the following assumptions. 

jff, jS, v, p, r, <j>, will all be zero identically (B IO) 

The equation for longitudinal motion, based on (B.5) and (B.6) are 



and the remaining Euler equation (B.4) 

6 = q 


(B.12) 


It is more convenient to select state variables a, V, q, 6, instead of u, w, q, 6 (B.4). 
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Using the relationship (Fig. B-4, /5 = 0) 

U = V cos a 
W = V sin a 
and V 2 = u 2 + w 2 

By taking derivatives with respect to time, we get 

uw-uw V a + sin a u 

a = =>• w = 

y2 COS a 


(B.13) 


Using the expression for w in (B.ll), after some algebraic manipulation, we get the equation for a: 

a T T x T_ 

a = q + — (cos 6 cos a + sin 6 sin a) - - sm a + cos a 

4 V mV mV mV 

In a similar fashion we derive expressions for V, q, and 6: 


11 iy * 2 

V = -g sin 6 cos a + g sin a cos 6 - — + — cos a + — sin a 

mm m 

. , qScCm , D ( ., iC0Stt , , jSina) , (< lS taa»«,co.o) * - «,,T 2 f I3) 

lyy lyy lyy lyy 

^ = q 

By substitution of the form of aerodynamic coefficients (B.9) and defining 


q n = -t z cos a + f x sina 
qi 2 = f z sina + f x cosa 


(B.14) 


and some more simple algebra, we get 


a = 


i Sc p 
1 +p — C 


_ p Sc 


4m 


_ SC L 


4m « 


(B.15) 


T T 

— — sin a + — — cos a 
mV mV 
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aS T x T z • 

V = -g sin 7 - -5_ Cr, + — cos a + — sin a 
m ° m m 


r. = P s 


q = 


21 


yy 


fC m< ,* qil C Do *q 12 C L JV 2 * ^Vq(SC mi _*q 12 Cg 


+ P Sc (cC ra . +q i2 J 


41 


yy 


ih.£15c, 

4m “ 


1 - C L ] JL cos 7 - 
4m L< i V 


P S C L T 

V - — — sin a 

2m mV 


+ cos a 

mV 


V + 


^ zt Tx “ ^xf T z 


yy 


* = q 


These are the four equations of motion used in this work. 


B.5 Values of Constants and Aerodynamic Coefficients 

The values of the various constants for this airplane were taken from [B.4], 


Constants: 


aircraft 

— McDonnel-Douglas F-18 Fighter 

mass 

m 

= 1035.308 slugs 

weight 

w 

= 33,3101b 

moments of inertia: 



Ixx 

= 23,000 slug-ft 2 


b. 

= 151,293 slug-ft 2 


Izz 

= 169,945 slug-ft 2 


bz 

= -2,971 slugs -ft 2 

wing area 

s 

= 400 ft 2 

MAC 

c 

= 11.52 ft 

wing span 

b 

= 37.42 ft 

geometry 

*x 

= -0.297 ft 


b 

= 0 



= 0.233 ft 


(B.15) 
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‘it 

= -19.37 ft 


ty( 

= 0 


t zt 

= 0.233 ft 

dynamic pressure 

q 

= 1/2 p V 2 

air density 

p 

= 0.0023709 o 


Aerodynamic Coefficients 

Reference [B.4] gives the aerodynamic coefficients for this aircraft in the form of curve fits to 
wind tunnel data. Values are given for the following range of parameters: 
angle of attack -10° < a < 90° 

side slip -20° < j8 < 20° 

mach number 0.2 < M < 2 

altitude 0 < h < 60,000 ft 

aileron deflection -25° < 6a < 25 

rudder deflection -30° < 5r < 30° 

stabilator deflection -24° < 6h < 10.5° 

throttle setting 30° (idle) < 6p < 131° (full afterburner) 

The coefficients are given in the form of piecewise arc tangent functions for discrete values of 
the parameters. Interim values are then interpolated. 

Example: 

C Mo - pitching moment in Eq. (B.9): 

C^(a,M = 0.6, 6h = 10.5°, h = 15,000 ft) = Cj^ X 6 (a) 

X 6 (a) = (.26/2.75) tan 1 (-(a°-5) 1/10) + (-.39/2.75) tan 1 ((a°-l) 1/8) 

+ (.8/2.75) tan* 1 ((a°-5) 1/13) + (.70/2.75) tan' 1 (-(a°-10) 1/65) 

+ (1.2/2.75) tan’ 1 (a°-49) 1/15) + (2.1/2.75) tan' 1 (-(-a°-69) 1/15) 

+ (-.45/2.75) tan’ 1 K<x°-77) 1/2) - 3.98 

Similar expressions are given for all the other relevant aerodynamic coefficients. This data is 
sufficient for the longitudinal case and also for small lateral angles. Ref. [B.4] also shows the accuracy 
of the curve filled data by comparing it to the actual wind tunnel results. In most cases the agreement 
between the curve filled and actual data is excellent. Small deviations can be expected in the high angle 
of attack range (J3 > w-70°). 
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The values of the aerodynamic coefficients are given for the following range of parameters: 

C D Sh = 10.5°, 0°, -5°, -24° 

h = 15 k' Sh = 10.5° M = 0.6, 0.9 

Sh = -24° M = 0.6, 0.9 

C Lq M = 0.6 h = 15 k' 

C L<i M = 0.6 h = 15 k' 

C mo M = 0.6, 0.9 Sh = 10.5°, 5°, 2°, 0°, -5°, -12.5°, -24° h = 15 k' 

C m M = 0.6 h = 15 k' 

C m . M = 0.6 h = 15 k' 

Further aerodynamic data can be found in Ref. [B.5] and [B.6] (manufacturers data) and several 
more reports which are in our hands. 

B.6 Future Developments 

The method of derivation for the equations of motion shown here can be used for the lateral mode 
also. In order to check the longitudinal controllers which were developed and to develop lateral 
controllers (if necessary) for combined pitch, yaw, and roll maneuvers which are needed in modern air 
combat, it will be necessary to develop a GDOF aerodynamic model, including GDOF equations of 
motion. This will be done gradually, by first constraining the lateral movements to small side slip angles, 
as was done by several references — for example, Safanov [B.7] discussing the Herbst [B.5] maneuver 
and Ostroff [B.9]. One can do this in two ways — either develop the full GDOF equations, from [B.5] 
and [B.6], and add a small /S constraint to the result, or use the small constraint from the beginning to 
get a somewhat simpler equation system, and move on to the full case later on. For the full GDOF case 
the data given in [B.4] must be widened by using the full data available in [B.5] and [B.6]. The group 
has already done preliminary work in this direction. A convenient method for showing results on real 
time simulations is given in Ref. [B.10]. This method uses spherical mapping to transform the various 
angular relations to a two-dimensional plane. Future lateral results may be displayed in this form. 
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